Gradient Formulation for the Stability of DC-Microgrids

TL;DR

This paper introduces a novel non-linear stability analysis method for dc-microgrids using gradient systems, enabling efficient equilibrium calculation and stability conditions through convex optimization, with practical validation via simulations.

Contribution

It presents a comprehensive stability analysis framework for dc-microgrids based on gradient system theory, addressing multiple aspects simultaneously.

Findings

Demonstrates existence and uniqueness of equilibrium

Provides numerical methods for equilibrium calculation

Establishes conditions for global stability

Abstract

This paper presents a non-linear stability analysis for dc-microgrids in both, interconnected mode and island operation with primary control. The proposed analysis is based on the fact that the dynamical model of the grid is a gradient system generated by a strongly convex function. The stability analysis is thus reduced to a series of convex optimization problems. The proposed method allows to: i) demonstrate the existence and uniqueness of the equilibrium ii) calculate this equilibrium numerically iii) give conditions for global stability using a Lyapunov function iv) estimate the attraction region. Previous works only address one of these aspects. Numeric calculations performed in cvx and simulations results in Matlab complement the analysis and demonstrate how to use this theoretical results in practical problems.

Tables

Table 1. TABLE I: Parameters of the constant power terminals and results of each model

Node	$k_i(\%)$	$p_i\left(kW\right)$	$c_i\left(\mu F\right)$	Model 1	Model 2
1	1	-80	690	0.9893	0.0958
2	1	-70	560	0.9851	0.0958
3	1	80	630	0.9917	0.0001
4	1	-50	500	0.9881	0.0710
5	1	-70	690	0.9954	0.0958
6	1	80	640	1.0047	0.0001
7	1	100	600	1.0129	0.0001
8	1	50	630	1.0026	0.0001
9	1	-90	690	0.9973	0.0818

Equations

$$c_i\frac{d{v}_i}{dt}=\frac{p_i+k_i(1-v_i)}{v_i}-g_{0i}{v}_0-\sum _{i=1}^n{g}_{ij}{v}_j$$

$$\frac{{\partial }^{2}W}{\partial x^2}⪰\mu I_N$$

$$M(x)\frac{dx}{dt}=-\frac{\partial W}{\partial x}$$

$$\frac{d\mathcal{V}}{dt}$$

$$\begin{array}{r}W(v)=-(p_i+k_i)\ln (v_i)+\sum _{i=1}^n{g}_{oi}{v}_0{v}_i+\\ \frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{g}_{ij}{v}_i{v}_j+\frac{1}{2}\sum _{i=1}^n{k}_i{v}_i^2\end{array}$$

$$\frac{{\partial }^{2}W}{\partial v^2}=(P+K)X+G+K$$

$$\Omega_{x}=\left\{\begin{array}[]{rl}(P+K)X+G+K&\succeq\mu I_{N}\\ X&=diag(x)\\ x_{i}\geq\;\tilde{x}_{i}&\geq 0\end{array}\right\}$$

$$\underset{v_i}{\min }W(v)$$

$$\underset{\alpha ,x_i}{\max }$$

$$(P+K)X+G+K$$

$$X$$

$$x_i$$

Full text

I Introduction

DCmicrogrids and dc distribution are promising technologies for integrating solar panels, batteries and fuel-cells among other components that operate in dc. Most of these components are integrated through a power electronic converter that controls a constant power resulting in a non-linear circuit [1]. In particular, constant power loads can introduce a negative resistance effect which in turns generates transient stability problems. Rigorous methodologies for stability analysis are required in this new context.

Several stability methodologies have been proposed for systems with ad-hoc controls [2]. However, the conventional droop control with a constant power model is the most common approach for control and stabilization of microgrids [3]. Most of the stability studies on this type of controls are based on linearization (i.e small signal stability). Transient stability in these generalized models is still a challenge due to their non-linear behavior; even finding the equilibrium point can be a challenge.

This paper proposes a methodology for transient stability of dc-microgrids based on the study of gradient systems. Although these type of systems have a rich and general theory, we are interested in a particular type, namely, gradient systems with a strongly convex function (we will call these as strongly-convex-gradient systems). We demonstrate the existence and uniqueness of the equilibrium as well as the conditions for global stability. In addition, we show a simple method for calculating this equilibrium and also estimating the region of attraction. Being a convex problem, we can guarantee convergence of the algorithms. Our model is simple enough to be tractable computationally but showing the main interaction between components of the nonlinear circuit. Just as in the case of the second order model for power systems applications. The proposed analysis shows a surprising connection between dynamical systems and convex analysis. This connection is explored from a practical point of view, since the stability analysis is transformed into a series of convex optimization problems that can be solved numerically. To the best of the author’s knowledge, there is not applications of the approach proposed in this paper.

The existence of the equilibrium was analyzed in [4] for one constant power load and generalized in [5] for several loads. A different approach was presented in [6] and [7] based on the convergence of the power flow. In the later methods, it was demonstrated that conventional algorithms such as Gauss and Newton’s methods converge to a unique equilibrium point under well defined conditions.

From the stability point of view, several studies have been presented for the small signal case [8][9][10]. However, transient stability studies are required in order to increase the accuracy of the study and consider the non-linear behavior of the grid. A recent review of transient stability analysis in microgrids can be found in [11].That review showed the necessity of systematic methods for large-scale stability analysis and reduced order models of the grid. Our method fulfills these conditions. In [12] a stabilization method was proposed for a dc-microgrid in which all the terminals were connected to the same bus-bar with only one equivalent constant power load. Droop controls were considered only on the sources. Our method consider the topology of the grid with different constant power loads and droop control in both the constant voltage and constant power terminals. In fact, constant power loads have been the main concern of recent stability analysis such as [13] and [14]. However, most of these studies are developed for ad-hoc controls requiring a detailed model of the converter. None of these approaches reveals the gradient characteristics of the model.

The use of convex analysis is also a contribution of this paper. Although this type of analysis is usually consider in linear matrix inequalities, that type of analysis is linear whereas the method presented here is nonlinear. In [14] a semidefinite programming methodology was proposed by formulating a Lure problem with quadratic bounds. That formulation allowed to estimate the region of attraction in grids with constant power loads. However the topology of the grid was limited to a unique bus-bar and the analysis is basically linear.

The rest of the paper is organized as follows: Section II presents the dynamical model of the dc-microgrid considering constant power terminals. Section III describes the stability analysis based on the use of gradient systems with strongly convex functions. Section IV explain how the stability problem is transformed into convex optimization models which allow to determine the attraction region, equilibrium point and under-voltage limit. Simulation results are presented in Section V followed by conclusions, appendix and references.

II Problem definition

Let us consider a dc-microgrid with droop control that is expected to operate whether grid-connected or in island-mode. The master node is represented by [math] and maintains a constant voltage $v_{0}$; the rest of the nodes are represented by $N=\left\{1,2,\cdots n\right\}$ and maintain a constant power which can be positive or negative. The grid is purely restive and is represented by the admittance matrix which includes linear loads eliminated by a kron’s reduction [15]. The model of each constant power terminal is depicted in Fig 1. It includes the capacitive effect of the converter and the droop control. This model is widely used in different applications including dc-microgrids [16] and multiterminal HVDC transmision [17].

The dynamics of the grid can be represented as follows:

[TABLE]

where $v_{i}$ are nodal voltages, $g_{ij}$ are entries of the admittance matrix and $c_{i}$ is the capacitance of the converter. Notice this is a non-linear dynamical system due to the presence of constant power devices. The model is the same for grid connected or island operation. The only difference is that under island operation, the switch is opened and hence $g_{i0}=0$.

On the other hand, each converter is equipped with a droop control with constant $k_{i}$. We assume that $c_{i}>0$ and $k_{i}\geq 0$ meaning that we can have terminals that do not participate in the primary control ($k_{i}$ can be zero). In addition, the graph that represents the grid is connected and hence $G=[g_{ij}]\succ 0$ (i.e it is positive definite). We allow transients generated by disconnection of the master node or abrupt changes in generation or demand. However, we assume the grid is not under short circuit and hence $v_{i}>0$ (our model would produce an infinite current for a shorcircuit in the capacitor).

It is important to notice that this model is general enough and allows to represent different components according to their functionality. For example, a solar panel is just a constant power device with $p_{i}>0$ while a constant power load is the same model but with $p_{i}<0$. As we will demonstrate below, constant power generation does not create big stability problems being the constant power loads the main source of instabilities. A simple justification of this fact, is that constant power loads generates a negative resistance effect that induce instability.

III Stability Analysis

III-A Convex analysis and gradient systems

We are interested in using convex optimization for analyzing stability in a particular type of dynamical systems, namely, gradient systems with strongly convex functions. Therefore, we will use some results from convex optimization restricting our analysis to continuously differentiable functions, in order to be used as Lyapunov candidates. Following this idea, let us start by this classical result

Theorem 1**.**

Consider a non-empty convex set $\Omega\in\mathbb{R}^{n}$ and a twice differentiable function $W:\Omega\rightarrow\mathbb{R}$ such that

[TABLE]

where $\mu$ is a real number such that $\mu>0$, $\partial^{2}W/\partial x^{2}$ is the Hessian matrix of $W$ and $I_{N}$ is the identity matrix of size $n$. We say that $W$ is a strongly convex; for this type of function there exists a global minimum $\tilde{x}\in\Omega$ and this minimum is unique.

Proof.

See [18] ∎

Remark 1**.**

There is a more general theorem which includes strictly convex functions (see [18] for more details). However, Theorem 1 is general enough for our dynamical problem and avoids the proliferation of unnecessary definitions and technicalities.

Remark 2**.**

It is well known that convex functions have a global optimum. However, for the case strongly convex functions, this optimum is not only global but also unique.

Let us connect this result with the dynamical system theory by proposing the following definition

Definition 1**.**

A strongly-convex-gradient system is a dynamical system that can be represented as

[TABLE]

where $M\succ 0$ (i.e positive definite) and $\partial H/\partial x$ represents the gradient of a strongly convex function $W:\Omega\in\mathbb{R}^{n}\rightarrow\mathbb{R}$.

With this simple definitions we can present our first result:

Theorem 2**.**

Consider a strongly-convex-gradient system described by (3). Then, there is a unique equilibrium point $\tilde{x}\in\Omega$ and this equilibrium is asymptotically stable. In addition $\Omega$ is an estimation of the atraction region.

Proof.

Consider the gradient system (3) with a strongly convex function $W$. The equilibrium point is given by the points in which $\partial W/\partial x=0$ which correspond to the global minimum of $W$. This minimum is also unique due to Theorem 1, guaranteeing the existence and uniqueness of the equilibrium point. For the stability analysis, consider a Lypunov function $\mathcal{V}(x)=W(x)-W(\tilde{x})$ which evidently fulfills the conditions for stability, namely $\mathcal{V}(\tilde{x})=0$, $\mathcal{V}(x)>0,\;\forall x\in\Omega-\left\{\tilde{x}\right\}$ and

[TABLE]

Asymptotic stability is directly demonstrated by invoking LaSalle’s invariant principle principle.

∎

This simple result allows to transform the problem of stability into a convex optimization problem. This is an advantage since we can calculate numerically the equilibrium point and estimate the region of attraction as will be demonstrated in the next sections.

III-B Dc-microgrids as gradient systems

Dc-microgrids can be represented as strongly-convex gradient systems. Let us formalize this by the following lemma that can be demonstrated by simple substitution and corresponding derivation:

Lemma 1**.**

The system given by (1) can be represented as a gradient system with gradient $W$ given by

[TABLE]

and $M=diag(c_{i})\succ 0$. In addition, the Hessian matrix of $W$ is given by

[TABLE]

where $P=diag(p_{i})$, $K=diag(k_{i})$ and $X=diag(1/v_{i}^{2})$

In order to demonstrate stability, it is enough to establish that exists a convex set $\Omega$ in which $W$ is strongly convex. This criteria is presented in the following theorem which constitute the main theoretical result of this paper:

Theorem 3**.**

A dc-microgrid with the conditions of Lemma 1 has an asymptotically stable equilibrium point if there exists a non-empty set containing the equilibrium $\tilde{x}$ which is given by

[TABLE]

with $\mu>0$. This set is an estimation of the region of attraction.

Proof.

It is easy to see from Lemma 1 that the function $W$ is strongly convex if the Hessian given by (6) is such that $\partial^{2}W/\partial v\succeq\mu I_{N}$. This is evidently true if the set (7) exists. This is an estimation of the region of attraction since it is an invariant set. ∎

Remark 3**.**

Recall that one of our main assumptions is that no capacitor is in short-circuit meaning that $v_{i}>0$. Therefore, the change of variables $x_{i}=1/v_{i}^{2}$ is a bijection, i.e we can obtain directly the value of $v_{i}$ from the value of $x_{i}$ and vice versa. We define two spaces $\Omega_{x}=\left\{x\in\mathbb{R}^{n}_{+}\right\}$ and $\Omega_{v}=\left\{v\in\mathbb{R}^{n}_{+}\right\}$ in order to simplify our nomenclature.

This theorem, allows a simple methodology for analyzing stability in dc-microgrids using convex optimization. For example, a direct result is obtained for the generation case:

Corollary 1**.**

For the generation case $(p_{i}>0)$ there is a unique equilibrium point and the region of attraction is the entire first quadrant (i.e $v_{i}>0$)

Proof.

Notice that for the generation case $p_{i}>0$ and hence we can increase $X$ as large as desired in order to fulfill the generalized inequality. The set of constraints is fulfilled even for $x_{i}\rightarrow\infty$ which correspond to the entire fisrt quadrant in both, the space $x$ and the space in $v$. ∎

This corollary only formalize a result that is well known for practical applications: generation case is not an stability issue being the constant power loads the real problem. In the next subsections, we analyze the general case with both constant power generation and constant power loads.

IV Stability analysis using convex optimization

Theorem 3 allows to define different convex optimization models for the stability analysis. In the following subsections we present different optimization models according to the desired analysis.

IV-A Equilibrium point

The existence and uniqueness of the equilibrium point on dc microgrids has been analyzed previously for different authors (see for example [6] and [7]). Here, we use Theorem 3 and define a simple optimization model

Model 1** (equilibrium point).**

[TABLE]

Hence, the equilibrium point of (1) can be obtained as $\tilde{v}=\text{argmin}\left\{W(v)\right\}$, given that $\Omega_{x}$ exists and contains the point $\tilde{x}_{i}=1/\tilde{v}_{i}^{2}$. The optimization problem can be solved by gradient or Newton’s methods guaranteeing uniqueness and convergence if $W$ is strongly-convex. The gradient method has a linear convergence while the newton method has a quadratic convergence (see [18] for more details).

IV-B Region of attraction

For the general case, the region of attraction can be estimated by maximizing the hypervolume of $\Omega_{x}$. However, the shape of $\Omega_{x}$ is a spectrahedron which can be quite complicated for being used in practical applications. For the sake of simplicity, we inscribe a hypercube $\Gamma_{x}$ inside $\Omega_{x}$ defined by $\Gamma_{x}=\left\{x\in\Omega_{x}:x_{i}\geq\alpha\right\}$. Then, our objective is to maximize the hypervolume of $\Gamma_{x}$ by maximizing $\alpha$. Therefore, an inner approximation of the region of attraction can be obtained by solving the following optimization model

Model 2** (region of attraction).**

[TABLE]

Notice this is a convex optimization problem. If the result is a positive $\alpha$ then the equilibrium is asymptotically stable and the values of $x_{i}$ define the limits of the hyper cube that represents an estimation of the region of attraction (in Section VI we show that this estimation is very accurate in practice).

It is possible to return to the original values by using the transformation from $\Omega_{x}$ to $\Omega_{v}$. In practice, this set can be considered as a limit for a voltage stability analysis and as an under-voltage protection, i.e the component must be disconnected for voltages lower than these voltages. Notice that Model 2 does not depend on the voltage $v_{0}$ and hence it can be used during voltage sags in the main grid (we will talk more about this issue in the results section).

V Results

The proposed methodology was evaluated in the 10-node dc-microgrid depicted in Fig 2. Nominal values of the grid are 380V/1kW and each line segment has a resistance of $1.5m\Omega/m$. The rest of the parameters are given in Table I.

Before considering the entire grid, led use study the case in which only load $p_{4}$ and generation $p_{7}$ are connected. The motivation of such a basic case, is that the dynamical model is in $\mathbb{R}^{2}$ allowing graphical representations (see [19] for a nullcline analysis on the plane for this simplified case).

The equilibrium point was obtained by Model 1 using the Newton’s method giving $v_{4}=0.8903$ and $v_{7}=0.9427$. Theorem 1 guarantee these values are unique in the convex set $\Omega$ which defines the region of attraction. This region can be estimated by Model 2 obtaining the minimum voltages $v_{min}=[0.182,\;0.002]$ which define an open subset in $\mathbb{R}^{2}$ as depicted in Figure 3; this figure depicts also the equilibrium point and three trajectories from different initial conditions calculated by the function ode45 of Matlab. As expected, initial points inside the estimation of the region of attraction go directly to the equilibrium. A fourth trajectory starting in the point $(0.1,0.1)$ (outside of the region of attraction) was calculated but not depicted in the figure since this trajectory was unstable and the function ode45 diverged.

Let us consider now the general case with 7 constant power terminals: 3 constant power generations (nodes 3,7,8) and 4 constant power loads (nodes 2,3,5,9). The main results of each model are given in Table I. Let us analyze each one of them.

The equilibrium point and the estimation of the region of attraction were obtained by Models 1 and 2 respectively. The estimation of the region of attraction is a hypercuve with the same interpretation as in the previous case. Now, this region can be used as criteria for voltage protection. In order to evaluate the under-voltage protection criteria, let us consider the case of a voltage sag in the main converter; the grid starts from an initial condition of $v_{0}=1pu$ but the voltage is reduced until $v_{0}=0.2pu$ for $t=0.05s$. As result, nodal voltages reach values bellow $v_{min}$ and the system becomes unstable (instability is diagnosed numerically since the function ode45 of Matlab diverge). This problem was avoided by disconnecting terminals with values below $v_{min}$. Figure 4 shows this type of protection for Load $p_{2}$ (similar behaviour is obtained for the rest of the terminals). Notice that the region of attraction was calculated considering the entire grid, but the criteria for the under-voltage protection is completely local, meaning that only its own measure of voltage is required.

These numerical simulations are available in [20] in order to reproduce the results and evaluate other possible configurations.

VI Conclusions

A convex based methodology for transient stability analysis of dc-micrids was presented. This methodology was based on a gradient representation of the grid and the general observation that the energy function is strongly-convex. The dynamical problem was then transformed into optimization problems in which Model 1 allowed to find the equilibrium point and Model 2 estimated the region of attraction. The later is also used in practice as a criteria for under-voltage protection. The main result was demonstrated mathematically using the Lyapunov method under well defined consideration regarded the convexity of the Lyapunov function. A complete set of simulations was performed in order demonstrate the practical applicability of the methodology.