Stability analysis of smart power grids with market dynamics in a port-

(1)

faculty of science and engineering

Stability analysis of smart power grids with market dynamics in a port-

Hamiltonian framework

Master Project Mathematics

December 2017

Student: J.W. Zijlstra

First supervisor: Prof. dr. A.J. van der Schaft Second supervisor: Prof. dr. C. De Persis

(2)

1 Introduction

1.1 Motivation and problem formulation

Provisioning energy for smart power grids

The interest in smart grids has increased steadily over the last decade. The traditional power grid consists of power plants that are responsible for gener- ating energy, which is transmitted and distributed to the end-users, i.e. the consumers. In this current day and age we are no longer dependent on one single energy source and shifted the focus nowadays on renewable energy. One of these renewable energies is solar power, where an increasing number of households install solar panels on their roof tops. As a consequence, households are no longer merely consumers but also producers of energy. Excess of energy can be sold and distributed to others who have a shortage. Hereto, households along with others (e.g. offices) have become active components in the sharing of power, which gives rise to smart power grids.

The provisioning of energy in these smart power grids has become increas- ingly difficult as generators operate more often at their limits and transmission lines can experience congestion. The use of real-time dynamic pricing as a control method can alleviate some of these problems [1]. In the regions where power supply cannot keep up with demand real-time dynamic pricing can incite consumers to change their usage. In this manner, consumers and producers fairly share their utilities and costs associated with the power consumption and production. The challenge lies in obtaining this in an optimal manner while the grid operates within its limits. This is called the social welfare problem [2], [3].

The market mechanics can be used to determine the optimal power dispatch.

Hereto, the market update process is dynamically linked to the physical response of the physical power network dynamics.

Extension

This subject has been broadly examined by literature. Stegink’s contribution consists of a rigorous and unifying passivity-based stability analysis using a port-Hamiltonian framework [5]. This framework is important as both the physical model as well as the dynamic pricing controllers admit a port- Hamiltonian representation and therefore are passive systems. Interconnection of port-Hamiltonian systems is power-preserving which implies passivity of the interconnected model.

For the physical model Stegink uses a third-order model to describe the synchronous generators, which is commonly referred to as the flux-decay model.

Stegink’s interconnection of the physical model with market dynamics gave rise to a couple of open questions which we will address later [5].

Conventionally, a lower-order model, i.e. the swing equations, is used in literature. As the same problems arise for the simpler dynamics and moreover the analysis is similar, we restrict ourselves to the lower-order model in trying to solve the open questions stated by Stegink. As such, the work proposed in this thesis forms an extension of the work done by [Stegink, 2017] which leads to the following problem formulation:

(5)

Thesis problem formulation

Stability analysis of smart grids with market dynamics in a port-Hamiltonian framework.

We work along the same lines as [Stegink, 2017]. One of the open problems that Stegink proposed is generalizing the type of graphs for which local asymptotic stability of the interconnected model with transmission costs and line congestion is proven.

The physical power flow cannot be directly controlled, however the interconnected model admits a virtual power flow on which one can place constraints.

Through the virtual flow we are able to control the physical flow if a one-to- one relation between the two is obtained. Stegink’s interconnected model with transmission costs and line congestion does not allow for this relation between flows to take place. As such, another aim of this thesis is to obtain power flow equality for all types of connected graphs.

1.2 Literature overview

As mentioned before, we mainly work along the same lines as [Stegink, 2017].

However, others have shown interesting extensions and approaches for the subject worth mentioning.

In [4] a fourth-order model for the physical model is used which includes tur- bine and exciter dynamics. An interconnection with a simple model describing the market dynamics is made.

There are several ways to solve an optimization problem. Commonly, one solves a general optimization problem like the social welfare problem by applying the primal-dual gradient method. The literature on the gradient method is quite vast [7], [8], [9]. Also in power grids this method is often applied to design distributed controllers, see for instance [10], [11], [12], [13], [14] or [15].

Stegink’s contribution shows that the real-time dynamic pricing model obtained when applying the gradient method can be represented as a port-Hamiltonian system [18], such that a port-Hamiltonian framework also arises for the market dynamics [5], [6]. Another method which is quite often used is the distributed averaging (integral) controller. This controller can be used for trying to solve optimal frequency regulation problems [21], [22].

Comparison of both controllers

Stegink discusses in [6] the pros and cons of both methods. A downside for using the distributed averaging (integral) controller is that it requires quadratic cost and utility functions whereas the primal-dual gradient method is able to deal with general convex utility functions.

A common feature of both controllers is that there exists freedom in choosing the communication graph, as long as the graph is connected.

In this thesis we will restrict ourselves to the primal-dual gradient controller.

(6)

Linear power flows have been extensively examined by literature, e.g. [10], [11], [12], [13]. However, [14], [15] show how optimal power dispatch can be achieved in power networks with non-linear power flows using gradient-method- based controllers. Downside of the controllers being used for non-linear power flows require the network topology to be a tree [14], [15].

1.3 Main contributions

Following in the foot steps of [Stegink, 2017] we interconnect the physical model with the market dynamics. However, as mentioned before a second-order model, i.e. the swing equations [6], is used to describe the physical model contrary to the third-order model Stegink adheres to [5]. As both interconnections give rise to the same open questions mentioned before, the shortcomings of the interconnected model are being examined and analyzed. The points where the interconnected model falls short are used as the building blocks for deriving a new interconnected model. This new model is also formed in a power-preserving port-Hamiltonian framework. Both algebraically and numerically we will investigate the properties of this newly formed model. Numerically we restrict ourselves to two network graphs. First, tendencies of the dynamics are examined for a simple 3-cycle, whereafter an extension to a bigger (random) graph is made.

Downside of the model is that non-convex constraints are used and as such the optimization problem is no longer convex. Hereto, the Karush-Kuhn-Tucker conditions are not sufficient to prove that the solution set is optimal [19], [20].

2 Preliminaries

2.1 Notation

The matrices occurring in this thesis consist of real elements, which for a given matrix A, is written as A ∈ R^n×m. Here n denotes the rows associated with the matrix, whereas m represents the columns. We use A > 0 (A ≥ 0) to indicate the matrix is positive (semi-)definite. In the case of vectors, v ≥ 0 denotes the element-wise inequality. The symbol 1 is used to represent a vector consisting of entries equal to 1. For an given vector η ∈ R^m we write sin(η) ∈ R^m to represent an element-wise sine function. Lastly, we write ∇f (x) and ∇²f (x) to represent respectively the gradient and Hessian of a twice-differentiable function f : Rⁿ→ Rⁿ.

Systems that admit the following structure:

˙

x = [J (x) − R(x)]∇H(x) + g(x)u

y = g(x)^T∇H(x) (1)

are called port-Hamiltonian input-state-output systems or port-Hamiltonian systems for short [18].

(7)

Here J (x) is some antisymmetric matrix, i.e. J (x) = −J (x)^T, and R(x) a positive semidefinite matrix, R(x) ≥ 0. The inputs of the system are represented by u ∈ R^m, the state variables by x ∈ Rⁿ, time-derivative ˙x ∈ Rⁿ and y ∈ R^m describes the output of the system.

The function H(x) is called the Hamiltonian and describes the energy associated to the system. Lastly, g(x) ∈ R^n×m is any given (matrix) function of state variable x which is linked to the input u of the system.

2.2 Convex optimization

The aim of an optimization problem is to find some x^∗∈ X such that it mini- mizes the problem, i.e. f (x^∗) = min{f (x) : x ∈ X }. Here X ⊂ Rⁿ denotes some feasible set and f (x) : Rⁿ → R is the objective function of the problem. The optimization problem becomes a convex optimization problem if X is a convex set and f (x) is a convex function defined on Rⁿ.

A constrained optimization problem of the form

minimize f (x) (2)

subject to gi(x) ≤ 0, i = 1, . . . , m (3) is called convex if the functions f, g1. . . gm: Rⁿ→ R are all convex functions.

One can also use the “standard form”, which consists of the following three parts:

• A convex objective function f (x) : Rⁿ→ R which will be minimized over the variable x,

• Inequality constraints of the form g_i(x) ≤ 0,where the functions g_i are convex and,

• Equality constraints of the form hi(x) = 0, where the functions hi are affine.

Affine functions are constraints of the form hi(x) = a^T_ix + bi, where ai∈ Rⁿ is a column-vector and bi∈ R a real number. Therefore, one can write a convex minimization (optimization) problem (CO) as

minx f (x) (4)

s. t. gi(x) ≤ 0, i = 1, . . . , m (5) hi(x) = 0, i = 1, . . . , p. (6) Remark: Note that every equality constraint hi(x) = 0 can be represented by a pair of inequality constraints hi(x) ≤ 0 and −hi(x) ≤ 0. For this reason, equality constraints seem to be redundant; however, it can be beneficial to treat them separately in practice. Observe that equality constraints have to be affine.

If hi(x) is convex, hi(x) ≤ 0 is convex, but −hi(x) ≤ 0 is concave. Therefore,

(8)

the only way for hi(x) = 0 to be convex is for hi(x) to be affine.

An important and sufficient condition for strong duality to hold for convex problems is the Slater condition, which is useful for checking optimal solutions of optimization problems [19], [20]. More on this in chapter 8.

Definition 2.1. A vector (point) x⁰∈ X is called a Slater point of the convex optimization problem (CO) (4) if

g_i(x⁰) < 0, for all i where g_i is nonlinear, g_i(x⁰) ≤ 0, for all i where g_i is linear.

If a Slater point exists, we say (CO) satisfies the Slater condition.

2.3 Lyapunov’s criterion in stability analysis

In this section we will state a method for checking stability of equilibrium points for a particular system consisting of ordinary differential equations (ODEs).

The method being considered is LaSalle’s invariance principle and the theory of Lyapunov, also referred to as the Lyapunov stability criterion. The criterion makes use of a Lyapunov function, V (x), which can be used to check stability of ODEs. The Lyapunov function has an analogy to the potential function in classical mechanics.

For a particular system ˙x = f (x), with equilibrium at x = ¯x, we consider a function V (x) : Rⁿ→ R where,

• Property 1: V (x) = 0 iff x = ¯x

• Property 2: V (x) > 0 for all x 6= ¯x

• Property 3: ˙V (x) = _dt^dV (x) ≤ 0 for all x 6= ¯x

If all three properties hold, the function is said to be a Lyapunov candidate for checking stability. If ˙V (x) is negative definite, the global asymptotic stability of ¯x is a consequence of Lyapunov’s criterion. LaSalle’s invariance principle gives a criterion for asymptotic stability in the case when ˙V (x) is only negative semidefinite. If V (x) > 0 and ˙V (x) ≤ 0 when x 6= ¯x hold only for x in some neighborhood Υ of the equilibrium ¯x, and the set { ˙V (x) = 0}T Υ does not contain any trajectories of the system besides the trajectory x(t) = ¯x, t ≥ 0, then the local version of LaSalle’s invariance principle states that ¯x is locally asymptotically stable.

(9)

3 Physical model and maximizing the social wel- fare

3.1 Power network model

Consider a power grid consisting of n buses. The network is represented by a connected and undirected graph G = (V, E ). Here V = {1, . . . , n} represents the nodes, i.e. the set of buses, and E ⊂ V × V are the edges interconnecting the nodes. These edges can be seen as the transmission lines connecting the buses.

The incidence matrix D of the network will be constructed by arbitrarily assigning a + and − to the ends of edge k. In this manner the incidence matrix of the resulting directed graph is given by,

Dik=







+1 if i is the positive end of k

−1 if i is the negative end of k

0 otherwise.

The following figure depicts an example of 3 buses interconnected to each other.

2 1 3

The incidence matrix corresponding to above network is given by,

D =





−1 0 −1

1 −1 0

0 1 1





Each bus represents a control area and is assumed to have a controllable power supply and demand. Therefore, each bus acts as a producer and consumer, i.e.

a prosumer. The dynamics at each bus is assumed to be given by,

δ˙_i= ω_i Miω˙i= − X

j∈N_i

ViVjBijsin(δi− δj) − Aiωi+ Pgi− Pdi (7)

which are referred to as the swing equations [6], [16], [22]. One can extend the dynamics by including the transient internal voltage of a synchronous generator.

These newly obtained dynamics is commonly known as the flux-decay model [5], [16], [17]. This thesis’ focus will lie on the swing equations to simplify simulations. However, forthcoming model analysis and calculations are still applicable to the flux-decay model when minor alterations are made [5].

(10)

An overview of all symbols and their meaning is listed in the next table.

δi Voltage angle ω_i^b Frequency

ωⁿ Nominal frequency

ωⁱ Frequency deviation ωi:= ω_i^b− ωⁿ Pdi Power demand

Pgi Power generation Mi Moment of inertia

Ni Set of buses connected to bus i A_i Asynchronous damping constant B_ii Self-susceptance

B_ij Negative of the susceptance of transmission line (i, j)

Throughout this thesis a number of assumptions will be made that are often made in power engineering [16].

Assumption 3.1.

• There is no conductance on the line, i.e the lines are lossless. In the case of high voltage lines connecting different control areas this assumption is valid.

• The network is operating around the nominal frequency ωⁿ. This implies that all our buses operate around the synchronous frequency (50 Hz).

• The nodal voltages Vi are constant.

On every edge in the network a voltage angle difference will occur. These differences between the buses will be denoted by η := D^Tδ. Furthermore, the angular momenta in the power grid will be given by p := M ω. Here ω represents the frequency deviations from the synchronous frequency, i.e ω = ω^b− 1ωⁿ, and M = diag_i∈V{M_i} is the diagonal matrix consisting of all moments of inertia in the power grid.

The swing equations (7) can be put in a more compact form. Let Γ = diag_k∈E = γk and γk = BijViVj, where k corresponds to the edge between nodes i and j. Hereto the swings equations can be compactly written as,

˙

η = D^Tω

M ˙ω = −DΓsin(η) − Aω + P_g− P_d (8) where D ∈ R^n×m is the incidence matrix corresponding to the network, A = diag_i∈VAi and column vectors Pg, Pd ∈ Rⁿ consisting respectively of elements Pgi and Pdi. Observe that Γsin(η) represents the power flow on the edges of the network, where sin(η) ∈ R^m.

(11)

The Hamiltonian Hp(η, p) associated to the swing equations is defined by Hp(η, p) = 1

2p^TM⁻¹p − 1^TΓcos(η) (9) The first term of the Hamiltonian represents the kinetic energy of the system, whereas 1^TΓcos(η) corresponds to the energy stored in the transmission lines, where cos(η) ∈ R^m. With help of the Hamiltonian, (8) can be put in a port- Hamiltonian form,

˙ xp = ˙η

˙ p

=

0 D^T

−D −A

∇Hp(η, p) +0 0 I −I

up

y =0 I 0 −I

∇Hp(η, p) = M⁻¹p

−M⁻¹p

= ω

−ω

(10)

The input upis defined as up= col(Pg, Pd). The above derived port-Hamiltonian system is hereafter referred to as the physical model.

3.2 Social welfare problem

A utility function U (Pd) can be assigned to the power consumption Pd and likewise a cost function C(Pg) to the power production Pg. We then define the social welfare by S(Pg, Pd) = U (Pd) − C(Pg). The objective is to maximize the social welfare while achieving zero frequency deviation. The assumption is made that C(Pg) is a strictly convex function and U (Pd) a strictly concave function. Therefore, the social welfare function is a strictly concave function in the variables Pg and Pd. In analyzing the equilibria of the compact swing equations (8) we will premultiply both left and right side by 1^T. One of the properties of an incidence matrix is that the column entries each add up to zero and as such, 1^TD = 0.

At equilibrium:

0 = − 1^TDΓsin(¯η) − 1^TA¯ω + 1^TP¯g− 1^TP¯d

= − 1^TA¯ω + 1^TP¯g− 1^TP¯d

Therefore, a necessary condition for zero frequency deviation, ¯ω = 0, is 1^TPd = 1^TPg, which is a mathematical way of stating that the total power generation must match the total demand at equilibrium.

From previous power balance one can state the following:

Theorem 3.2. Let D_c ∈ R^n×m^c denote an incidence matrix of some connected communication graph with m_c edges and n nodes. Then, a solution (P_g, P_d) of the power balance, 1^TP_d = 1^TP_g, exists if and only if there exists a v ∈ R^m^c satisfying D_cv − P_g+ P_d= 0.

(12)

Proof. (⇒) In the case of power balance we have 1^TPd = 1^TPg, and thus Pg− Pd ∈ ker 1^T ⇐⇒ Pg− Pd ∈ (span{1})^⊥. As Dc is a directed graph, each column sums up to 0. Thus, D_c^T sums each row up to 0 and therefore D^T_cx = 0 =⇒ x1 = x2 = ... = xn = α with α being some scalar (xi are the elements of the vector x). Therefore, ker D_c^T = span 1, which in turn implies im D^T_c = (span{1})^⊥. Combining earlier result we obtain P_g− Pd ∈ im Dc. Thus, ∃v ∈ im D_c such that D_cv = P_g− Pd.

(⇐) Assume D_cv − P_g + P_d = 0 holds true. Premultiplying latter equation by 1^T and observing 1^TD_cv = 0 one immediately obtains the desired result:

1^TP_d= 1^TP_g.

We will use latter results to show that the social welfare problem admits a natural convex optimization form, which we will address hereafter. From (4), (5), (6) maximizing the social welfare is equivalent to the following minimization problem:

min

(P_g,P_d,v)− S(Pg, Pd) = C(Pg) − U (Pd) s.t. − Dcv + Pg− Pd= 0,

(11)

which has a constrained optimization form, and due to the constraints being affine, the social welfare problem can be interpreted as a convex optimization problem [19], [20].

The Lagrangian associated to maximizing the social welfare (11) is given by,

L = C(Pg) − U (Pd) + λ^T(Dcv − Pg+ Pd), (12) where λ ∈ Rⁿdenotes the Lagrange multiplier. From the Lagrangian one can determine its corresponding Karush-Kuhn-Tucker (KKT)-conditions, which arise by taking the gradient of the Lagrangian and setting it to zero, i.e ∇L = 0

∇C( ¯P_g) − ¯λ = 0 (13a)

−∇U ( ¯Pd) + ¯λ = 0 (13b)

D_c^Tλ = 0¯ (13c)

Dcv − ¯¯ Pg+ ¯Pd= 0 (13d) Due to the optimization problem being convex and the Slater condition being satisfied, it follows that ( ¯Pg, ¯Pd, ¯v) is an optimal solution to (11) if and only if there exists λ such that (13) holds, a result from literature as strong duality holds [19], [20]. Chapter 8 explains duality in more detail. As we will see in upcoming section, the KKT-conditions give rise to the sought-after market dynamics.

(13)

4 Basic primal-dual gradient controller

The primal-dual gradient method is widely used in power engineering literature.

The method will act as a basis for constructing a dynamic pricing algorithm to represent the market dynamics, which will be used as a starting point for the controller. Applying the primal-dual gradient method [7], [8], [9] to the earlier derived KKT-conditions (13), one obtains:

τgP˙g= −∇C(Pg) + λ + u^g_c (14a) τ_dP˙_d= ∇U (P_d) − λ + u^d_c (14b)

τv˙v = −D^T_cλ (14c)

τλ˙λ = Dcv − Pg+ Pd (14d)

Here τc = blockdiag (τg, τd, τv, τλ) > 0 are parameters for designing the controller and determine the convergence speed of the dynamics. A “large” τc

results in a slow convergence, whereas a “small” τc speeds up the convergence rate. Additional inputs have been introduced uc = col(u^g_c, u^d_c) which will be specified later on. For the communication graph Dc that appears in (14c) &

(14d) one has freedom in choosing its structure. One can decide to take the same structure of the physical network, but one may also prefer for instance to have all-to-all communication for which the graph is complete.

The equations (14) have an economic interpretation. For the power generation dynamics (14a) one can observe that the producers try to minimize the cost function towards an equilibrium where the gradient of the cost function equals λ + u^g_c. Differently stated, each power producer aims to maximize his own profits, which occurs whenever their marginal cost is equal to the local price λi+ u^g_ci. The same analogy can be made for (14b) where this time it is the consumers looking to maximize their own utility but are amerced by the local price λi− u^d_ci. The dynamic (14c) describes the distribution of the pricing mechanism in the physical network, whereas v, obtained from (14d), represents a virtual power flow on the edges of the communication graph described by Dc. The word virtual arises as the flow described by v may differ from the physical power flow due to using a different (communication) graph. As mentioned before, uc = col(u^g_c, u^d_c) are additional inputs to the system. They describe penalties or prices that are assigned to the power producers and consumers respectively. The inputs are chosen in such a way that they compensate for the frequency deviation occurring in the physical network.

In the next section we will interconnect the markets dynamics (14) to the physical system to obtain a closed-loop model and values for u_c = col(u^g_c, u^d_c) are specified.

(14)

4.1 Interconnecting physical power network with market dynamics

In the last section it became clear that we are heading towards an interconnection between the physical network and accompanied market dynamics. To this end, we will start off this section with the derivation of this closed-loop model. Later on, the equilibria of this newly formed model will be addressed and we end this section comparing the differences that occur when either cyclic or acyclic graphs are being considered.

4.1.1 Derivation of closed-loop model

Recall the market dynamics from (14). As we are primarily interested in the system’s behaviour around steady-state, for ease, we will consider the case where the controller design parameters are given by τc = blockdiag(τg, τd, τv, τλ) = I, where I represents the identity matrix. For a complete model containing these parameters we refer to [5].

The dynamics (14) are a system of ordinary differential equations (ODEs).

Together with the physical model they can be represented as:

˙

η = D^Tω

M ˙ω = −DΓsin(η) − Aω + Pg− Pd

P˙g= −∇C(Pg) + λ + u^g_c P˙d= ∇U (Pd) − λ + u^d_c

˙v = −D^T_cλ

˙λ = Dcv − Pg+ Pd

(15)

However, defining the variables xc= (xg, xd, xv, xλ) a port-Hamiltonian representation for the market dynamics can be obtained, namely

˙ xc =







0 0 0 I

0 0 0 −I

0 0 0 −D^T_c

−I I D_c 0







∇Hc(xc) + ∇S(xc) +





 I 0 0 I 0 0 0 0





 uc

y_c=I 0 0 0

0 I 0 0

∇Hc(x_c) =P_g P_d

,

(16)

where uc= col(u^g_c, u^d_c) is the to be determined additional input and S(xc) is the social welfare problem stated in (11), as such

∇S(xc) =







∇P_gS(xc)

∇P_dS(xc) 0 0





 .

(15)

Strictly speaking, latter obtained form is an incremental port-Hamiltonian system as S is concave and therefore satisfies the passivity property:

(x1− x2)^T(∇S(x1) − ∇S(x2)) ≤ 0, ∀x1, x2∈ R^3n+m^c. (17) As the swing equations also admit a port-Hamiltonian form one can interconnect both systems by an appropriate choice for the input functions u_c and u_p. Let uc = −yp = −ω

ω

and up = yc = Pg

Pd

, then the following interconnected port-Hamiltonian form is obtained:

˙ x =







0 D^T 0 0 0 0

−D −A I −I 0 0

0 −I 0 0 0 I

0 I 0 0 0 −I

0 0 0 0 0 −D_c^T

0 0 −I I D_c 0







∇H(x) + ∇S(x), (18)

where x is the combination of x_p= col(η, M ω) = col(η, p) and x_c = col(P_g, P_d, v, λ) and H(x) represents the total Hamiltonian of the interconnected system, given by H(x) = H_p(x_p) + H_c(x_c). Here H_c(x_c) =¹₂x^T_cx_c, a natural quadratic choice for our controller Hamiltonian and H_p the swing equations Hamiltonian as defined earlier, H_p(η, p) = ¹₂p^TM⁻¹p − 1^TΓcos(η).

The system (18) also lends itself a compact port-Hamiltonian form, namely

˙

x = (J − R)∇H(x) + ∇S(x) (19)

The matrix J − R has the structure

J − R =







0 D^T 0 0 0 0

−D −A I −I 0 0

0 −I 0 0 0 I

0 I 0 0 0 −I

0 0 0 0 0 −D_c^T

0 0 −I I Dc 0





 ,

where J = −J^T and R = R^T ≥ 0.

In the following section we will use Lyapunov analysis as a method to prove (local) stability of the closed-loop system.

(16)

4.2 Lyapunov stability analysis of the closed-loop model

To investigate the stability of the previously obtained model (18), we have to examine its equilibria. Hereto, we define the following equilibrium set for the closed-loop model:

X1= {¯x|¯x is an equilibrium of (18)} (20) The elements of this set are the desired equilibria as they satisfy by definition the KKT-conditions (13) and due to (13d) zero frequency is achieved for the physical model (10).

For the physical model we will assume the following for equilibrium points:

Assumption 4.1. There exists an equilibrium ¯x ∈ X₁ that satisfies η ∈ (−π/2, π/2)^mand η ∈ im D^T. Moreover, the π/2−constraint implies that the Hessian of the Hamiltonian evaluated at this equilibrium point is positive definite, i.e. ∇²H(¯x) > 0.

The assumption η ∈ (−π/2, π/2)^m is made as this acts as a security constraint for power grid stability [22].

Next theorem shows a neighborhood exists for which convergence to the optimal points, i.e. equilibria takes place. The theorem states additionally that the convergence is to a point. We omit this part of the proof as it is quite cumbersome, however interested readers are referred to [Stegink, 2017] for the complete proof.

Theorem 4.2. For every ¯x ∈ X₁satisfying Assumption 4.1 there exist a neighborhood Υ around ¯x where all trajectories x satisfying (18) with initial conditions in Υ converge to the set X₁.

Proof. We will use Lyapunov’s criterion for checking stability of the equilibrium points of the closed-loop system (18). A natural Lyapunov candidate to consider for the closed-loop system is the Hamiltonian. Recall from section 2.3 that in order for a function V (x) to be a Lyapunov function it has to satisfy all three properties. The first property required V (x) = 0 iff x = ¯x. However, the Hamiltonian does not necessarily have to be 0 at x = ¯x. As such, we make use of the shifted Hamiltonian around ¯x to ensure the first property is being satisfied.

The shifted Hamiltonian around ¯x is defined as,

H(x) = H(x) − (x − ¯¯ x)^T∇H(¯x) − H(¯x)

Note that the shifted Hamiltonian is not globally bounded as the terms (x − ¯x)^T∇H(¯x) − H(¯x) lead to a linear contribution which is not bounded.

As such, the shifted Hamiltonian only allows local stability analysis and thus requires LaSalle’s invariance principle.

Hereto, consider the Lyapunov function V (x) = ¯H(x)

(17)

Due to the definition of the shifted Hamiltonian property 1 is automatically met (locally). The second property, V (x) > 0 for all x 6= ¯x in a small neighbourhood around ¯x, comes from local strict convexity around ¯x, because of Assumption 4.1 we have ∇²V (x) = ∇²H(x) = ∇¯ ²H(x) > 0. In order to see that the third property also holds, one can observe that

V (x) =˙ d

dtV (x) = ∇ ¯H(x)^Tx˙

For the latter property to hold true we need the latter equation to be less or equal to zero for all x 6= ¯x in Υ. Rewriting ˙x for the shifted Hamiltonian version of equation (19), one obtains using ∇H(x) = ∇ ¯H(x) + ∇H(¯x) and 0 = (J − R)∇H(¯x) + ∇S(¯x)

˙

x = (J − R)∇H(x) + ∇S(x)

= (J − R)(∇ ¯H(x) + ∇H(¯x)) + ∇S(x)

= (J − R)∇ ¯H(x) + ∇S(x) − ∇S(¯x) Hence the time-derivative ˙V admits to,

V (x) = ∇ ¯˙ H(x)^Tx˙

= ∇ ¯H(x)^T((J − R)∇ ¯H(x) + ∇S(x) − ∇S(¯x))

For the first term ∇ ¯H(x)^T(J − R)∇ ¯H(x), observe that this is the same as

−∇ ¯H(x)^TR∇ ¯H(x) by making use of the anti-symmetric property J = −J^T. The second term one can notice that

∇ ¯H(x) = ∇H(x) − ∇H(¯x) =







M⁻¹(p − ¯p) Γsin(η) − Γsin(¯η)

Pg− ¯Pg

Pd− ¯Pd

v − ¯v λ − ¯λ







=





M⁻¹(p − ¯p) Γsin(η) − Γsin(¯η)

xc− ¯xc





Moreover, recalling that C(P_g) and U (P_d) are respectively strictly convex and concave we use the structure of

∇S(x) =





 0 0

∇P_gS(x)

∇P_dS(x) 0 0





 ,

together with the form of ∇ ¯H(x) to obtain the following time-derivative of V (x) V (x) = −ω˙ ^TAω + (x − ¯x)^T(∇S(x) − ∇S(¯x)) ≤ 0. (21)

(18)

Bearing in mind that A > 0 and that S(x) is strictly concave in Pg and Pd, we have that both terms in the time-derivative of V (x) are non-positive.

Moreover, equality in (21) holds iff Pg= ¯Pg, Pd= ¯Pdand ω = 0.

Due to ˙V (x) being non-positive there exist a compact sublevel set Υ of V (x) where the initial conditions in Υ converges to the largest invariant S contained in Υ ∩ {x|Pg= ¯Pg, Pd= ¯Pd, ω = 0}.

On such invariant set the dynamics of the (shifted) closed loop model (15)- (18) changes into







˙ η

˙ p P˙g

P˙d

˙v

˙λ







=







˙ η 0 0 0

˙v

˙λ







=







0 D^T 0 0 0 0

−D −A I −I 0 0

0 −I 0 0 0 I

0 I 0 0 0 −I

0 0 0 0 0 −D_c^T

0 0 −I I D_c 0













Γsin(η) − Γsin(¯η) M⁻¹(p − ¯p)

P_g− ¯P_g P_d− ¯P_d v − ¯v λ − ¯λ







+







0 0

−∇C(Pg) + ∇C( ¯Pg)

∇U (Pd) + ∇U ( ¯Pd) 0

0







(22) This gives rise to the following system of equations:

1. ˙η = D^TM⁻¹(p − ¯p) = 0 ⇐⇒ η = constant 2. λ − ¯λ = 0 ⇐⇒ λ = ¯λ

3. ˙λ = Dc(v − ¯v)

4. ˙v = −D^T_c(λ − ¯λ) = 0, by λ = ¯λ → v = constant

Thus, λ = ¯λ and η, v are constant. Hence, x initialized in Υ converges to S ⊂ X1 as t → ∞. Therefore, each point in the largest invariant set where V = 0 corresponds to an equilibrium, which concludes the proof.˙

Variations in the controller design

There are several variations and extensions one can make to the social welfare problem (11). For instance, one can include nodal power constraints on the power consumption and production. Stegink explains how this can be incorpo- rated into the original problem [5], [6]. However, this thesis’ focus lies on another extension one can make to the optimization problem (11), where transmission costs and line congestion are included. In the next section we will introduce this modified social welfare problem and derive the accompanied closed-loop system.

(19)

4.3 Social welfare problem including line congestion and transmission costs

Define the modified social welfare problem as S(P_g, P_d, v) = U (P_d) − C(P_g) − C_T(v), where C_T(v) is a convex function representing the power transmission cost. Moreover, we can bound the virtual power flow on each transmission line as a security constraint to prevent overloads.

As such, the following modified optimization problem (11) is obtained

min

(P_g,P_d,v)− S(Pg, Pd, v) = C(Pg) + CT(v) − U (Pd) s.t. Dv − P_g+ P_d= 0

− κ ≤ v ≤ κ

(23)

Here κ ∈ R^m satisfies the element-wise inequality κ > 0. Note that we have taken the communication graph Dc to be identical to the topology of the physical network D. As a result, κ denotes the security constraint placed on the (virtual) power flow along each transmission line. As we have done for the original social welfare problem (11) a system of ODEs can be derived from (23).

As before, the starting point is the Lagrangian of the system, given by

L = C(P_g) + C_T(v) − U (P_d) + λ^T(Dv − P_g+ P_d) + µ^T₊(v − κ) + µ^T₋(−κ − v)

with Lagrange multiplier λ ∈ Rⁿ and µ+, µ₋ ∈ R^m_≥0. From the Lagrangian, the corresponding KKT optimality conditions are obtained

∇C( ¯Pg) − ¯λ = 0 (24a)

− ∇U ( ¯P_d) + ¯λ = 0 (24b)

∇CT(¯v) + D^Tλ + ¯¯ µ+− ¯µ₋ = 0 (24c)

D¯v − ¯P_g+ ¯P_d= 0 (24d)

− κ ≤ ¯v ≤ κ (24e)

¯

µ₊, ¯µ₋ ≥ 0, ¯µ₊^T(¯v − κ) = 0, ¯µ^T₋(−κ − ¯v) = 0 (24f) Since our constraints are linear, we have that the modified problem satisfies Slater condition. As such, since the problem (23) is convex, ( ¯P_g, ¯P_d, ¯v) is an optimal solution of (23) iff there exists λ ∈ Rⁿ and µ₊, µ₋ ∈ R^m_≥0 satisfying previous KKT-conditions.

If we apply the gradient method in a similar way as for the original KKT- conditions we obtain the following closed-loop system of ODEs:

(20)

˙

η = D^Tω (25a)

M ˙ω = −DΓsin(η) − Aω + P_g− Pd (25b)

P˙g= −∇C(Pg) + λ − ω (25c)

P˙d= ∇U (Pd) − λ + ω (25d)

˙v = −∇C_T(v) − D^Tλ − µ₊+ µ (25e)

˙λ = Dv − Pg+ Pd (25f)

˙

µ+= (v − κ)⁺_µ₊ (25g)

˙

µ−= (−κ − v)⁺_µ₋ (25h)

where,

(v − κ)⁺_µ

+=

v − κ if µ+> 0 max{0, v − κ} if µ+= 0 (−κ − v)⁺_µ₋ =

−κ − v if µ₋> 0

max{0, −κ − v} if µ₋= 0 i ∈ E

The equations (25a)-(25f) lend itself a (partial) port-Hamiltonian form, namely

˙

x = (J − R)∇H(x) + ∇S(x) + N µ N =0 0 0 0 −I 0

0 0 0 0 I 0

^T , (26)

with J, R and x the same as before (19). In the next section we will investigate the equilibrium conditions of (26). As for the original basic model which did not include line congestion or transmission costs a Lyapunov analysis of the equilibrium points will be made.

4.4 Equilibrium analysis for model including line conges- tion and transmission costs

For the closed-loop system stated in (25), the following equations (25b, 25f) are obtained at equilibrium

0 = −DΓsin(¯η) − A¯ω + ¯Pg− ¯Pd

0 = D¯v − ¯P_g+ ¯P_d,

Using the power balance 1^TP¯g = 1^TP¯d for zero frequency deviation the latter equations simplify to:

D¯v = DΓsin(¯η) (27)

(21)

Next we will distinguish two scenarios such that we can say more about latter obtained result.

4.4.1 Scenario 1: Acyclic graphs

If the graph of D is taken such that it resembles a spanning tree (acyclic) then kerD = {0}. As such, taking a spanning tree as the topology for D results in the following equilibrium state: Γsin(¯η) = ¯v. Hence, one can directly link the controller variable v to the physical power flow of the network at steady state.

4.4.2 Scenario 2: Cyclic graphs

Upon investigating the general case where cycles are included in our graphs, the kernel of the incidence matrix D does not trivially equal 0 anymore. However, we are interested if some of the previous acyclic results where we could directly link v to the physical power flow still hold true.

In order to show the impact of cycles in our graphs, we will distinguish two cases. One for the linear approximation of the closed-loop system and one for the complete model. Upon investigating both cases a simple 3-cycle graph will be considered:

2 1 3

which has the following corresponding incidence matrix:

D =





−1 0 −1

1 −1 0

0 1 1



 Case 1: Cyclic graphs in linear approach

Recall from Assumption 3.1 that all our buses operate around the synchronous frequency. Typically in power engineering it is assumed that the voltage angle differences are small which allows for a linearization study of the closed-loop system, sin(η) = sin(D^Tδ) ≈ D^Tδ. We will assume the power flow v admits a similar form as the physical flow, i.e. Γsin(η) ≈ ΓD^Tδ. Hereto, we introduce a virtual voltage angle ξ, giving rise to the virtual power flow form v = ΓD^Tξ.

As such, investigating steady state conditions of (27), will lead to the following result:

−DΓD^Tδ + DΓD^Tξ = 0, which can be rewritten as:

DΓD^T(δ − ξ) = 0 ⇔ δ − ξ ∈ ker DΓD^T

(22)

With help of the next theorem we will show that latter result is equivalent to ΓD^Tδ = ΓD^Tξ, and thus Γsin(η) ≈ ΓD^Tδ = v for the linear case.

Theorem 4.3. The following is true for any incidence matrix D of a connected and directed graph:

ker DΓD^T = ker D^T = span 1

Proof. In order to prove this equality, we will show that ker DΓD^T ⊃ ker D^T and ker DΓD^T ⊂ ker D^T. Moreover, both statements turn out to be equal to span 1.

(⊂) For any x ∈ ker D^T, it follows that D^Tx = 0 ⇒ DΓD^Tx = 0. Thus, x ∈ ker DΓD^T and hence ker D^T ⊂ ker DΓD^T.

(⊃) For any x ∈ ker DΓD^T, we have DΓD^Tx = 0. Upon multiplying both sides with x^T, we obtain x^TDΓD^Tx = 0. Taking D^Tx = y one can simplify this to get y^TΓy = 0. Using the positive definite property of Γ, the previous statement holds only true when y = 0 ⇒ D^Tx = 0. As such, x ∈ ker D^T, which leads to ker DΓD^T ⊂ ker D^T.

Combining both results, proves our desired statement ker DΓD^T = ker D^T. As D is a directed graph, each column sums up to 0. Thus, D^T sums each row up to 0 and therefore D^Tx = 0 ⇔ x1= x2= ... = xn= α with α being some scalar (xiare the elements of the vector x). Therefore, ker D^T = span 1, which in turn leads to the desired outcome δ − ξ ∈ ker DΓD^T = span 1 ⇔ D^Tδ = D^Tξ.

Algebraic 3-cycle example

We will quickly demonstrate above mentioned result, D^Tδ = D^Tξ algebraically and use it for the complete model. As such, we are interested whether ker DΓD^T equals ker D^T = span 1.

DΓD^T =





−1 0 −1

1 −1 0

0 1 1









γ1 0 0 0 γ2 0 0 0 γ3









−1 1 0

0 −1 1

−1 0 1





=





γ1+ γ3 −γ1 −γ3

−γ1 γ1+ γ2 −γ2

−γ3 −γ2 γ2+ γ3



 For each x ∈ ker DΓD^T we have the following:





γ1+ γ3 −γ1 −γ3

−γ1 γ1+ γ2 −γ2

−γ3 −γ2 γ2+ γ3



x =



 0 0 0



, and by using Gaussian elimination one obtains:





−γ1 γ1+ γ2 −γ2

−γ3 −γ2 γ2+ γ3

0 0 0







 x1

x2

x3



=



 0 0 0



 (28)

(23)

Solving for x1gives:

x₁= γ₂x₃− (γ1+ γ₂)x₂

−γ1

One can easily check that substituting x1 in (28) leads to x2= x3 and thus:

x₁=γ₂x₂− (γ1+ γ₂)x₂

−γ1

= −γ1x₂

−γ1

= x₂= x₃ Hence ker DΓD^T = span 1.

Case 2: Cyclic graphs in complete model

For the linear case we could show that the virtual voltage angle differences should be equal to its physical counterpart, i.e. D^Tξ = D^Tδ. Once again we will let v = ΓD^Tξ. Using this for the equilibrium condition (27) results in,

DΓD^Tξ = DΓsin(¯¯ η)

DΓ D^Tξ − sin(D¯ ^Tδ) = 0¯ (29) The main question that arises is whether the following is true:

D^Tξ − sin(D¯ ^T¯δ) = 0,

as linearization leads to D^Tξ = D¯ ^Tδ, the desired result. Like we did for the¯ linear case we will use the 3-cycle example to see whether earlier stated question is true.

Algebraic 3-cycle example

From the particular 3-cycle incidence matrix we obtain,

D^Tξ =¯



 ξ¯2− ¯ξ1

ξ¯3− ¯ξ2

ξ¯3− ¯ξ1



, and

sin(D^T¯δ) =





sin(¯δ₂− ¯δ₁) sin(¯δ₃− ¯δ₂) sin(¯δ₃− ¯δ₁)



 So do we have the following:



 ξ¯2− ¯ξ1

ξ¯3− ¯ξ2

ξ¯3− ¯ξ1





=?





sin(¯δ2− ¯δ1) sin(¯δ3− ¯δ2) sin(¯δ3− ¯δ1)



 (30)

Generally speaking latter equality does not hold true. In section 5 we will address this in more detail numerically.

In the next section we will perform a Lyapunov analysis on the extended model for acyclic networks, as a one-to-one relation between physical and virtual power flow was shown for this particular case.

(24)

4.5 Lyapunov analysis on the extended closed-loop model for acyclic networks

The Lyapunov analysis provided in this section will only be applicable to acyclic networks as we are guaranteed a direct link between the controller variable v to the physical power flow, i.e. ¯v = Γsin(¯η).

Like we did for the closed-loop model, to investigate the stability of the extended model (25), we have to examine its equilibria. Hereto, we define the following equilibrium set for the extended closed-loop model:

X2= {(¯x, ¯µ)|(¯x, ¯µ) is an (isolated) equilibrium of (25)} (31) Theorem 4.4. Let the network topology be acyclic and let (¯x, ¯µ) ∈ X2 be an (isolated) equilibrium satisfying Assumption 4.1. Then all trajectories (x, µ) of (25) initialized in a small neighborhood Υ around (¯x, ¯µ) converge asymptotically to (¯x, ¯µ).

Proof. Like before, Lyapunov’s criterion will be used. Recall the shifted Hamil- tonian around ¯x,

H(x) = H(x) − (x − ¯¯ x)^T∇H(¯x) − H(¯x), which one can use to rewrite (26) as

˙

x = (J − R)∇ ¯H(x) + ∇S(x) − ∇ ¯S(x) + N ˜µ (32) where ˜µ = µ − ¯µ. Consider the Lyapunov function

V (x, µ) = ¯H(x) +1

2µ˜^T₊µ˜++1 2µ˜^T₋µ˜₋

Note that the first 2 Lyapunov properties are satisfied by using Assumption 4.1 and notice the similarity with previous used Lyapunov function, V (x) = ¯H(x).

In order to check whether

V =˙ d dtV is less or equal to zero, observe the following:

˜

µ^T₊(v − κ)⁺_µ₊≤ ˜µ^T₊(v − κ) = ˜µ^T₊(¯v − κ + ˜v) ≤ ˜µ^T₊v˜

˜

µ^T₋(−κ − v)⁺_µ₋≤ ˜µ^T₋(−κ − v) = ˜µ^T₋(−κ − ¯v − ˜v) ≤ −˜µ^T₋˜v (33) Using latter result, the time-derivative ˙V amounts to,

(25)

V =∇ ¯˙ H(x)^T((J − R)∇ ¯H(x) + ∇S(x) − ∇S(¯x)) − ˜µ^T₊˜v + ˜µ^T₋˜v + ˜µ^T₊(v − κ)⁺_µ

++ ˜µ^T₋(−κ − v)⁺_µ₋

= − ω^TAω + (x − ¯x)^T(∇S(x) − ∇S(¯x)) − ˜µ^T₊˜v + ˜µ^T₋˜v + ˜µ^T₊(v − κ)⁺_µ₊+ ˜µ^T₋(−κ − v)⁺_µ₋

≤ − ω^TAω + (x − ¯x)^T(∇S(x) − ∇S(¯x))

≤0

where equality holds iff ω = 0, Pg= ¯Pg, Pd= ¯Pd.

Due to ˙V being non-positive there exist a compact sublevel set Υ of V (x, µ) where the initial conditions in Υ converge to the largest invariant S contained in Υ ∩ {(x, µ)|P_g= ¯P_g, P_d= ¯P_d, ω = 0}.

On such invariant set our dynamics of the (shifted) closed loop model (25a)- (25f) changes into







˙ η

˙ p P˙g

P˙d

˙v

˙λ







=







˙ η 0 0 0

˙v

˙λ







=







0 D^T 0 0 0 0

−D −A I −I 0 0

0 −I 0 0 0 I

0 I 0 0 0 −I

0 0 0 0 0 −D^T

0 0 −I I D 0













Γsin(η) − Γsin(¯η) M⁻¹(p − ¯p)

P_g− ¯P_g P_d− ¯P_d v − ¯v λ − ¯λ







+







0 0

−∇C(Pg) + ∇C( ¯Pg)

∇U (Pd) + ∇U ( ¯Pd)

−∇CT(v) + ∇CT(¯v) 0





 +







0 0

−I I

0 0







µ+− ¯µ+

µ−− ¯µ−

(34)

This gives rise to the following system of equations:

1. ˙η = D^TM⁻¹(p − ¯p) = 0 ⇐⇒ η = constant 2. λ − ¯λ = 0 ⇐⇒ λ = ¯λ

3. ˙λ = D(v − ¯v), by ker D = 0 we have that v = ¯v

4. ˙v = −D^T(λ − ¯λ) − ∇CT(v) + ∇CT(¯v) − (µ+− ¯µ+) + (µ−− ¯µ−) = 0, by using λ = ¯λ and v = ¯v → µ+= ¯µ+, µ− = ¯µ−

Thus, λ = ¯λ, µ = ¯µ and η, v are constant. Hence, (x, µ) initialized in Υ converges to S ⊂ X2 as t → ∞. Therefore, each point in the largest invariant set where ˙V = 0 corresponds to an equilibrium, which concludes the proof.

Note that above mentioned proof is only valid for acyclic network topologies.

As such, in the next section we will look at the numerical results from both the basic and extended closed-loop model where acyclic graphs are included.

(26)

5 Numerical results closed-loop model

In this section we will look at the numerical results for both the basic unconstrained closed-loop model (15) as well as the constrained model (25) where transmission costs and line congestion have been taken into account.

For the numerical analysis we assume the following for the models:

• C(Pg) = ¹₂P_g^TQgPg

• U (P_d) = −¹₂P_d^TQ_dP_d+ b^TP_d

• CT(D^Tξ) =

(0 if unconstrained

1

2(D^Tξ)^TQT(D^Tξ) if constrained

The parameters for an arbitrary incidence matrix D ∈ R^n×m are given by:

Q_d = I, Q_T = 2I Qg = diag(1, 2, . . . , n)

A = 2I, Γ = I, b = col (1, 1.25, 1.5, 1, . . . , 1)

Let us first look at the numerical results of a single 3-cycle for both the constrained and unconstrained model.

5.1 Simple 3-cyclic graph

We will start off with the simplest case, i.e. a single cycle for the 3-node case,

2 1 3

with corresponding incidence matrix:

D =





−1 0 −1

1 −1 0

0 1 1





Firstly, we will examine the evolution for all state variables occurring in the system of ODEs (15) and compare them to the state variables of the extended model (25). The transmission cost parameter, QT, is kept constant, whereas the line congestion parameter κ will vary. Three cases will be investigated, namely κ = 0.1, κ = 0.2 and κ = 1 for fixed QT = 2I. The case κ = 1 is similar to the case where only transmission costs are considered, thus the power flow constraints are not active. The value for κ is large enough such that the security

Stability analysis of smart power grids with market dynamics in a port-