University of Groningen
Hybrid Interconnection of Iterative Bidding and Power Network Dynamics for Frequency
Regulation and Optimal Dispatch
Stegink, Tjerk; Cherukuri, Ashish; De Persis, Claudio; Van Der Schaft, Arjan; Cortes, Jorge
Published in:IEEE Transactions on Control of Network Systems DOI:
10.1109/TCNS.2018.2856404
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.
Document Version
Final author's version (accepted by publisher, after peer review)
Publication date: 2019
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Stegink, T., Cherukuri, A., De Persis, C., Van Der Schaft, A., & Cortes, J. (2019). Hybrid Interconnection of Iterative Bidding and Power Network Dynamics for Frequency Regulation and Optimal Dispatch. IEEE Transactions on Control of Network Systems, 6(2), 572-585. [8411329].
https://doi.org/10.1109/TCNS.2018.2856404
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.
Hybrid interconnection of iterative bidding and power network dynamics
for frequency regulation and optimal dispatch
Tjerk Stegink, Ashish Cherukuri, Claudio De Persis, Arjan van der Schaft, and Jorge Cort´es
Abstract—This paper considers a real-time electricity market involving an independent system operator (ISO) and a group of strategic generators. The ISO operates a market where genera-tors bid prices at which they are willing to provide power. The ISO makes power generation assignments with the goal of solving the economic dispatch problem and regulating the network frequency. We propose a multi-rate hybrid algorithm for bidding and market clearing that combines the discrete nature of iterative bidding with the continuous nature of the frequency evolution in the power network. We establish sufficient upper bounds on the inter-event times that guarantee that the proposed algorithm asymptotically converges to an equilibrium corresponding to an efficient Nash equilibrium and zero frequency deviation. Our technical analysis builds on the characterization of the robustness properties of the continuous-time version of the bidding update process interconnected with the power network dynamics via the identification of a novel LISS-Lyapunov function. Simulations on the IEEE 14-bus system illustrate our results.
I. INTRODUCTION
The dispatch of power generation in the grid has been traditionally done in a hierarchical fashion. Broadly speaking, cost efficiency is ensured via market clearing at the upper layers and frequency regulation is achieved via primary and secondary controllers at the bottom layers. Research on im-proving the performance of these layers has mostly developed independently from each other, motivated by their separation in time-scales. The increasing penetration of renewables poses significant challenges to this model of operation because of its intermittent and uncertain nature, see e.g., [2], [3]. At the same time, the penetration of renewables also presents an opportunity to rethink the architecture and its hierarchical separation towards the goal of improving efficiency and adap-tiveness. A key aspect to achieve the integration of different layers is the characterization of the robustness properties of the mechanisms used at each layer, since variables at the upper layers cannot be assumed in steady state any more at the lower ones. These considerations motivate our work on iterative bidding schemes combined with continuous physical network dynamics and the correctness analysis of the resulting multi-rate hybrid interconnected system.
A preliminary version of this work appeared as [1] at the American Control Conference.
This work is supported by the NWO Uncertainty Reduction in Smart Energy Systems (URSES) program, the ARPA-e Network Optimized Distributed Energy Systems (NODES)program and AFOSR Award FA9550-10-1-0499.
T. W. Stegink, C. De Persis and A. J. van der Schaft are with the Jan C. Willems Center for Systems and Control, University of Groningen, 9747 AG Groningen, the Netherlands.{t.w.stegink, c.de.persis, a.j.van.der.schaft}@rug.nl
A. Cherukuri is with the Automatic Control Laboratory, ETH Z¨urich.
cashish@control.ee.ethz.ch
J. Cort´es is with the Department of Mechanical and Aerospace Engineering, University of California, San Diego.cortes@ucsd.edu
Literature review: The integration of economic dispatch and frequency regulation in power networks has attracted increasing attention in the last decades. Many recent works [4], [5], [6], [7], [8], [9], [10], [11] envision merging the design of primary, secondary, and tertiary control layers for several models of the power network/micro-grid dynamics with the aim of bridging the gap between long-term optimization and real-time frequency control. In scenarios where generators are price-takers, the literature has also explored the use of market mechanisms to determine the optimal allocation of power gen-eration and to stabilize the frequency with real-time (locational marginal) pricing, see [12], [13], [14], [15]. Our present work shares with [16], [17] the use of dynamic iterative bidding schemes by the ISO, although in these works the setting is stochastic, the agents react in a price-taking manner, and their dynamics is assumed to be decoupled from one another. Instead, [18] proposes iterative bidding schemes where the generators are strategic, leading to efficient Nash equilibria where power generation levels minimize the total cost as intended by the ISO. Inspired by [18] and [15], our work [19] has shown that the integration with the frequency dynamics of the network can also be achieved in scenarios where generators are price-setters. However, the integration in [19] relies on a continuous-time model for the bidding process, where the frequency coming from the power network dynamics enters as a feedback signal in the negotiation process. Instead, we account here for the necessarily discrete nature of the bidding process and explore the design of provably correct multi-rate hybrid implementations that realize this integration. This involves the identification of a strictly decreasing Lyapunov function (compared to the weak Lyapunov function employed in [19]) that allows us to go beyond the characterization of asymptotic stability of the continuous-time dynamics and establish strong robustness guarantees. These guarantees set the basis for our design here of a provably correct, multi-rate hybrid algorithm implementation.
Statement of contributions: We consider an electrical power network consisting of an ISO and a group of strategic gener-ators. The ISO seeks to ensure that the generation meets the load with the minimum operation cost and the grid frequency is regulated to its nominal value. Each generator seeks to maximize its individual profit and does not share its cost function with anyone. The ISO operates the market, where generators bid prices at which there are willing to provide power, and makes power generation assignments based on the bids and the local frequency measurements. Our goal is to design mechanisms that ensure the stability of the intercon-nection between the ISO-generator bidding process and the physical network dynamics while accounting for the different nature (iterative in the first case, evolving in continuous time in
the second) of each process. Our starting point is a continuous-time bid update scheme coupled with the physical dynamics of the power network whose equilibrium corresponds to an effi-cient Nash equilibrium and zero frequency deviation. Our first contribution is the characterization of the robustness properties of these dynamics against additive disturbances. To achieve this, we identify a novel local Lyapunov function that includes the energy function of the closed-loop system. The availability of this function not only leads us to establish local exponential convergence to the desired equilibrium, but also allows us rigorously establish its local input-to-state stability properties. Building on these results, our second contribution develops a time-triggered hybrid implementation that combines the discrete nature of iterative bidding with the continuous nature of the frequency evolution in the power network. In our design, we introduce two iteration loops, one (faster) inner-loop for the bidding process that incorporates at each step the frequency measurements, and one (slower) outer-loop for the market clearing and the updates in the power generation levels, that are sent to the continuous-time power network dynamics. We refer to this multi-rate hybrid implementation as time-triggered because we do not necessarily prescribe the time schedules to be periodic. To analyze its convergence properties, we regard the time-triggered implementation as an approximation of the continuous-time dynamics and invoke the robustness properties of the latter, interpreting as a disturbance their mismatch. This allows us to derive explicit upper bounds on the length between consecutive triggering times that guarantee that the time-triggered implementation remains asymptotically convergent. The computation of these upper bounds does not require knowledge of the efficient Nash equilibrium. Simula-tions on the IEEE 14-bus power network illustrate our results. Notation: Let R, R≥0, R>0, Z≥0, Z≥1 be the set of real,
nonnegative real, positive real, nonnegative integer, and posi-tive integer numbers, respecposi-tively. For m ∈ Z≥1, we use the
shorthand notation [m] = {1, . . . , m}. For A ∈ Rm×n, we let
kAk denote the induced 2-norm. Given v ∈ Rn, A = AT ∈
Rn×n, we denote kvk2A:= vTAv. The notation 1 ∈ Rnis used
for the vector whose elements are equal to 1. The Hessian of a twice-differentiable function f : Rn→ R is denoted by ∇2f .
II. POWER NETWORK FREQUENCY DYNAMICS
Here we present the model of the physical power network that describes the evolution of the grid frequency. The network is represented by a connected, undirected graph G = (V, E ), where nodes V = [n] represent buses and edges E ⊂ V ×V are the transmission lines connecting the buses. Let m denote the number of edges, arbitrarily labeled with a unique identifier in [m]. The ends of each edge are also arbitrary labeled with ‘+’ and ‘-’, so that we can associate to the graph the incidence matrix D ∈ Rn×m given by Dik=
+1 if i is the positive end of edge k, −1 if i is the negative end of edge k, 0 otherwise.
(1)
A spanning tree T is as a connected acyclic subgraph of G that contains all vertices of G. The incidence matrix DT associated
to T is constructed as in (1). Each bus represents a control area and is assumed to have one generator and one load. Following [20], the dynamics at the buses is described by the following swing equations
˙δ = ω,
M ˙ω = −DΓ sin(DTδ) − Aω + Pg− Pd.
(2)
Here, Γ = diag{γ1, . . . , γm} ∈ Rm×m, γk = BijViVj,
where k ∈ [m] corresponds to the edge between nodes i and j. For the asymptotic stability analysis carried out in Sections IV and V, we assume that the load Pd is constant.
Table I specifies the meaning of the symbols used in the power
δ ∈ Rn (vector of) voltage phase angles
ω ∈ Rn frequency deviation w.r.t. the nominal frequency Vi∈ R>0 voltage magnitude at bus i
Pd∈ Rn constant power load
Pg∈ Rn power generation
M ∈ Rn×n≥0 diagonal matrix of moments of inertia
A ∈ Rn×n≥0 diagonal matrix of asynchronous damping constants
Bij∈ R≥0 negative of the susceptance of transmission line (i, j)
Table I: Parameters and state variables of model (2).
network model (2). The validity of this model relies on the following assumptions, which are standard in the literature on power network dynamics [20], [21], and we state here for completeness.
Assumption II.1. (Swing equation model): For the power network dynamics described by (2), the following hold:
• Lines are lossless, i.e., the conductance of all lines is zero. This is generally valid for transmission lines;
• Nodal voltages Vi are constant; • Reactive power flows are ignored;
• Network is balanced such that the three-phase network
can be analyzed by a single phase.
To avoid issues in the stability analysis of (2) due to the rotational invariance of δ, see e.g., [22], we introduce a new set of variables. To this end, consider an arbitrary spanning tree T of G and let DT be its associated incidence matrix.
Consider ϕ = DTTδ ∈ Rn−1 representing the voltage phase angle differences along the edges of this spanning tree. The physical energy stored in the transmission lines is given by
U (ϕ) = −1TΓ cos(DTD†TT ϕ), (3) where DT† = (DTTDT)−1DTT denotes the Moore-Penrose
in-verse of DT. By noting that DTD†TD = (I −1n11
T)D = D,
the system (2) in the (ϕ, ω)-coordinates takes the form
˙
ϕ = DTTω,
M ˙ω = −DT∇U (ϕ) − Aω + Pg− Pd.
(4)
In the sequel we assume that, for the power generation Pg=
¯
Pg, there exists an equilibrium col( ¯ϕ, ¯ω) of (4) that satisfies
DTD†TT ϕ ∈ (−¯ π2,π2)m. The latter assumption is standard and often referred to as the security constraint [20].
III. PROBLEM STATEMENT
In this section we formulate the problem statement and then discuss the paper objectives. We start from the power network model introduced in Section II and then explain the game-theoretic model describing the interaction between the ISO and the generators following the exposition of [23], [18].
The cost incurred by generator i ∈ [n] in producing Pgi
units of power is given by
Ci(Pgi) :=12qiPgi2 + ciPgi, (5)
where qi> 0 and ci≥ 0. The total network cost is then
C(Pg) := X i∈[n] Ci(Pgi) =12P T g QPg+ cTPg, (6)
with Q = diag{q1, . . . , qn} and c = col(c1, . . . , cn). Given
the cost (6) and the constant power loads Pd, the ISO seeks
to solve the economic dispatch problem minimize
Pg
C(Pg), (7a)
subject to 1TPg= 1TPd, (7b)
and, at the same time, regulate the network frequency to its nominal value. Since the function C is strongly convex, there exists a unique optimizer Pg∗ of (7). However, we assume that the generators are strategic and they do not reveal their cost functions to anyone, including the ISO. Consequently, the ISO is unable to determine the optimizer of (7). Instead, it determines the power dispatch according to a market clearing procedure in which each generator submits bids to the ISO.
We consider price-based bidding: each generator i ∈ [n] submits the price per unit electricity bi ∈ R at which it
is willing to provide power. Based on these bids, the ISO finds the power generation allocation that minimizes the total generator payment while meeting the load. More precisely, given the bid b = col(b1, . . . , bn), the ISO solves
minimize
Pg
bTPg, (8a)
subject to 1TPg= 1TPd. (8b)
The optimization problem (8) is linear and may in general have multiple (unbounded) solutions. Among these solutions, let Pgopt(b) = col(Pg1opt(b), . . . , P
opt
gn(b)) be the optimizer of
(8) the ISO selects given bids b. Knowing this process, each generator i aims to bid a quantity bi to maximize its payoff
Πi(bi, P opt gi (b)) := P opt gi (b)bi− Ci(P opt gi (b)). (9)
For an unbounded optimizer we have Πi(bi, ±∞) = −∞.
To analyze the clearing of the market, we resort to tools from game theory [24]. To this end, we define the inelastic electricity market game:
• Players: the set of generators [n].
• Action: for each player i ∈ [n], the bid bi∈ R. • Payoff: for each player i ∈ [n], the payoff Πi in (9).
We refer to the game as inelastic, as the load is not affected by the bids bi. For the bid vector we interchangeably use the
notation b ∈ Rn and (bi, b−i) ∈ Rn, where b−irepresents the
bids of all players except i. A bid profile b∗ ∈ Rn is a Nash
equilibrium if there exists an optimizer Pgopt(b∗) of (8) such
that ∀i ∈ [n], Πi(bi, P opt gi (bi, b∗−i)) ≤ Πi(b∗i, P opt gi (b ∗))
for all bi 6= b∗i and all optimizers P opt
gi (bi, b∗−i) of (8). In
par-ticular, we are interested in bid profiles that can be associated to economic dispatch. More specifically, a bid b∗ ∈ Rn is
efficient if there exists an optimizer Pg∗ of (7) which is also an optimizer of (8) given bids b = b∗ and
Pgi∗ = arg max
Pgi
{Pgib∗i − Ci(Pgi)} for all i ∈ [n]. (10)
A bid b∗ is an efficient Nash equilibrium if it is both efficient and a Nash equilibrium. At the efficient Nash equilibrium, the optimal generation allocation determined by (7) coincides with the production that the generators are willing to provide, maximizing their profit (9). Following the same arguments as in the proof of [23, Lemma 3.2], one can establish the existence and uniqueness of the efficient Nash equilibrium. Proposition III.1. (Existence and uniqueness of efficient Nash equilibrium): Let (Pg∗, λ∗) be a primal-dual optimizer of (7), then b∗ = ∇C(Pg∗) = 1λ∗ is the unique efficient Nash
equilibrium of the inelastic electricity market game.
In the scenario described above, neither the ISO nor the in-dividual strategic generators are able to determine the efficient Nash equilibrium beforehand. Our goal is then to design an online bidding algorithm where ISO and generators iteratively exchange information about the bids and the generation quan-tities before the market is cleared and dispatch commands are sent. The algorithm should be truly implementable, meaning that it should account for the discrete nature of the bidding process, and at the same time ensure that network frequency, governed by the continuous-time power system dynamics, is regulated to its nominal value. The combination of these two aspects leads us to adopt a hybrid implementation strategy to tackle the problem.
IV. ROBUSTNESS OF THE CONTINUOUS-TIME BID AND POWER-SETPOINT UPDATE SCHEME
In this section, we introduce a continuous-time dynamics that prescribes a policy for bid updates paired with the frequency dynamics of the power network whose equilib-rium corresponds to an efficient Nash equilibequilib-rium and zero frequency deviation. In this scheme, generators update their bids in a decentralized fashion based on the power generation quantities received by the ISO, while the ISO changes the generation quantities depending on both the generator bids and the network frequency. This design is a simplified version of the one proposed in our previous work [19]. The main contribution of our treatment here is the identification of a novel Lyapunov function that, beyond helping establish local exponential convergence, allows us to characterize the input-to-state stability properties of the dynamics. We build on this characterization later to develop our time-triggered hybrid implementation that solves the problem outlined in Section III.
A. Bidding process coupled with physical network dynamics Recall from Section III that given bid bi, generator i ∈
[n] wants to produce the amount of power that maximizes its individual profit, given by
Pgides:= arg max Pgi
{biPgi− Ci(Pgi)} = qi−1(bi− ci). (11)
Hence, if the ISO wants generator i to produce more power than its desired quantity, that is Pgi > Pgides, generator i will
increase its bid, and vice versa. Bearing this rationale in mind, the generators update their bids according to
Tb˙b = Pg− Q−1b + Q−1c. (12a)
Here Tb∈ Rn×n is a diagonal positive definite matrix. Next,
we provide an update law for the ISO depending on the bid b ∈ Rn and the local frequency of the power network. The ISO updates its actions according to
TgP˙g= 1λ − b + ρ11T(Pd− Pg) − σ2ω,
τλ˙λ = 1T(Pd− Pg).
(12b)
with parameters ρ, σ, τλ ∈ R>0 and where Tg ∈ Rn×n is a
diagonal positive definite gain matrix.
The intuition behind the dynamics (12b) is explained as follows. If generator i bids higher than the Lagrange multiplier λ (sometimes referred to as the shadow price [25]) associated to (8b), then the power generation (setpoint) of node i is decreased, and vice versa. By adding the term with ρ > 0, one can enhance the convergence rate of (12b), see e.g., [26]. We add the feedback signal −σ2ω to compensate for the
frequency deviations in the physical system. Interestingly, albeit we do not pursue this here, the dynamics (12) could also be implemented in a distributed way without the involvement of a central regulating authority like the ISO.
For the remainder of the paper, we assume that there exists an equilibrium ¯x = col( ¯ϕ, ¯ω, ¯b, ¯Pg, ¯λ) of (4)-(12) such that
DTD†T
T ϕ ∈ (−π/2, π/2)¯ m (cf. Section II). Note that this
equilibrium satisfies ¯ λ = 1 T(P d+ Q−1c) 1TQ−11 > 0, ω = 0,¯ ¯b = 1¯λ, ¯ Pg= Q−11¯λ − Q−1c, 1TP¯g = 1TPd. (13)
In particular, at the steady state, the frequency deviation is zero, the power balance1TP¯g= 1TPd is satisfied, and1¯λ =
¯
b = ∇C( ¯Pg), implying that ¯Pgis a primal optimizer of (7) and
¯
b is an efficient Nash equilibrium by Proposition III.1. Hence, at steady state the generators do not have any incentive to deviate from the equilibrium bid.
B. Local input-to-state (LISS) stability
While the ISO dynamics (12b) is a saddle-point dynamics of the linear optimization problem (8) (and hence, poten-tially unstable), we show next that the interconnection of the physical power network dynamics (4) with the bidding process (12) is locally exponentially stable and, furthermore, robust to additive disturbances. First, inspired by our previous work [19], we define the function
V (x) = U (ϕ) − (ϕ − ¯ϕ)T∇U ( ¯ϕ) − U ( ¯ϕ) +1 2ω TM ω + 1 2σ2(kb − ¯bk 2 Tb+ kPg− ¯Pgk 2 Tg + kλ − ¯λk 2 τλ). (14) where x = col(ϕ, ω, b, Pg, λ). Then the closed-loop system
obtained by combining (4) and (12) is compactly written as ˙ x = F (x) = Q−1AQ−1∇V (x) (15) with Q = QT = blockdiag(I, M,Tb σ, Tg σ, τλ σ) > 0 and A = 0 DT T 0 0 0 −DT −A 0 σI 0 0 0 −Q−1 I 0 0 −σI −I −ρ11T 1 0 0 0 −1T 0 .
By exploiting the structure of the system, we obtain the dissipation inequality
˙
V = 12(∇V (x))TQ−1(A + AT)Q−1∇V (x) ≤ 0. (16)
However, since R := −12(A + AT) is only positive
semi-definite, V is not strictly decreasing along the trajectories of (15). Nevertheless, one can employ this function, cf. [19], and invoke the LaSalle Invariance Principle to characterize the local asymptotic convergence properties of the dynamics. Here, we show that, in fact, the dynamics is locally input-to-state (LISS) stable, as defined in [27], and therefore robust to additive disturbances. Our key tool to establish this is the identification of a LISS-Lyapunov function, which in general is far from trivial for dynamics that involve saddle-point dynamics. To this end, consider the system
˙
x = F (x) + Bd (17)
for some B ∈ R4n×q and disturbance signal d ∈ Rq. While the function V defined in (14) does not qualify for being a LISS Lyapunov function, the next result shows that adding suitable cross-terms to it yields an LISS-Lyapunov function for the system (17).
Theorem IV.1. (LISS-Lyapunov function for the intercon-nected dynamics): Consider the interconintercon-nected dynamics (17) and define the function
W(x) = V (x) + 01(ϕ − ¯ϕ)TDT†M ω (18) −02 σ2 (b − ¯b) TT g(Pg− ¯Pg) −0σ23(λ − ¯λ)1 TT g(Pg− ¯Pg),
with parameters = col(0, 1, 2, 3) ∈ R4>0 andV given by
(14). Given the equilibrium ¯x = col( ¯ϕ, ¯ω, ¯b, ¯Pg, ¯λ) of (15), let
¯
η = DTD†T
T ϕ. For γ such that k¯¯ ηk∞ < γ < π2, define the
closed convex set
Ω = {x = col(ϕ, ω, b, Pg, λ) | DTDT†Tϕ ∈ [−γ, γ]m}. (19)
Then there exist sufficiently small such that W is an
LISS-Lyapunov function of (17) on Ω. In particular, there exist constants α, χ, c1, c2 > 0 such that for all x ∈ Ω and all
d satisfying kdk ≤ χkx − ¯xk, 1 2c1kx − ¯xk 2≤ W (x) ≤ 12c2kx − ¯xk2, (20a) (∇W(x))T(F (x) + Bd) ≤ −αkx − ¯xk2. (20b)
We refer to Appendix A for the proof of Theorem IV.1. Using the characterization (20) and [28, Theorem 4.10], each
trajectory of (15) initialized in a compact level set contained in Ω exponentially converges to the equilibrium ¯x corresponding to economic dispatch and the efficient Nash equilibrium. Moreover, we exploit the local ISS property of (17) guaranteed by Theorem IV.1 next to develop a time-triggered hybrid implementation.
V. TIME-TRIGGERED IMPLEMENTATION:ITERATIVE BID UPDATE AND MARKET CLEARING
In realistic implementations, the bidding process between the ISO and the generators is not performed continuously. Given the availability of digital communications, it is reason-able to instead model it as an iterative process. Building on the continuous-time bidding dynamics proposed in Section IV, here we develop a time-triggered hybrid implementation that combines the discrete nature of bidding with the continuous nature of the frequency evolution in the power network. We consider two time-scales, one (faster) for the bidding process that incorporates at each step the frequency measurements, and another one (slower) for the market clearing and updates of the power generation levels that are sent to the power network dynamics. We refer to this implementation as time-triggered because we do not necessarily prescribe the time schedules to be periodic in order to guarantee that the asymptotic stability properties are retained by the hybrid implementation.
A. Algorithm description
We start with an informal description of the iterative update scheme between the ISO and the generators, and the intercon-nection with the dynamics of the power network.
[Informal description]: The algorithm has two time indices, k to label the iterations on the bidding process and l to label the iteration in the market clearing process that updates the power setpoints. At each iteration l ∈ Z≥0, ISO and generators
are involved in an iterative process where, at each subiteration k, generators send a bid to the ISO. Once the ISO has obtained the bids and the network frequency measurements at time tl
k, it computes
the new potential generation allocations, denoted Pk+1
g ∈ Rn, and sends the corresponding one to
each generator. At the (k + 1)-th subiteration, each generator adjust its bid based on their previous bid and the generation allocation received from the ISO at time tlk+1. Once k = Nl∈ Z≥1 at time tlNl, the market is cleared, meaning that the bidding process is reset (i.e., k = 0), the power generations in the swing equations are updated according to the current setpoints PNl
g , and the index l moves to l + 1.
Figure 1 shows the two iteration layers in the update scheme. The evolution of the frequency occurs in continuous time according to (4). To relate iteration numbers with time instances on R, we consider time sequences of the form {{tl k} Nl k=0} ∞ l=0 which satisfy tlk− tl k−1> 0, t l+1 0 = t l Nl, ∀l ∈ Z≥0, ∀k ∈ [Nl]. (21)
Here Nl ∈ Z≥1 is the number of bid iterates before the
bidding process of market clearing instance l ∈ Z≥0stops, see
Figure 1. Algorithm 1 formally describes the iterative updates of the bidding process between the generators and the ISO.
Algorithm 1: ITERATIVE BID UPDATE AND MARKET
CLEARINGALGORITHM
Executed by: generators i ∈ [n] and ISO Data : time sequence {{tl
k} Nl
k=0} ∞
l=0; cost function
(5) for each generator i; frequency deviation ω(tl
k) at each time tlk and load Pd for ISO
Initialize : each generator i selects arbitrarily b0 i ≥ ci,
sets k = 0, l = 0, and jumps to step 6; ISO selects arbitrary Pgi0 > 0, λ0i > 0, sets k = 0, l = 0 and waits for step 8
1 while l ≥ 0 do
2 while k ≥ 0, k < Nldo
3 /* For each generator i: */ 4 Receive Pgik from ISO at tlk; Set
5 bk+1i = bki + (tlk+1− tlk)Tbi−1(Pgik − q−1i (bki + ci)) 6 Send bk+1i to the ISO; set k = k + 1
7 /* For ISO: */
8 Receive bki, ωi(tlk) from each i ∈ [n] at tlk 9 Set Pgik+1= Pgik + (tlk+1− tlk)Tgi−1(λk− bki −
σ2ωi(tlk) + ρ
P
i∈[n](Pdi− Pgik)) for all i ∈ [n]
λk+1= λk+tlk+1−t l k τλ P i∈[n](Pdi− Pgik) 10 Send Pgik+1to each i ∈ [n], set k = k + 1
11 end 12 Set Pgi(t) = PgiNl in (4) ∀i ∈ [n], ∀t ∈ [tlN l, t l+1 Nl+1) 13 Set b0i = b Nl i , P 0 gi = P Nl gi , λ 0 i = λ Nl i for each i ∈ [n] 14 Set l = l + 1, k = 0 15 end
For analysis purposes, we find it convenient to represent the dynamics resulting from the combination of Algorithm 1 and the network dynamics (4) as the time-triggered continuous-time system ˙ ϕ(t) = DTTω(t), M ˙ω(t) = −DT∇U (ϕ(t)) − Aω(t) + Pg(tl0) − Pd, Tb˙b(t) = Pg(tlk) − Q −1b(tl k) − Q −1c, (22) TgP˙g(t) = 1λ(tlk) − b(t l k) − σ 2ω(tl k) + ρ11 T(P d− Pg(tlk)), τλ˙λ(t) = 1T(Pd− Pg(tlk)), for t ∈ [tlk, tlk+1) ⊂ [t0l, tl+10 ), l ∈ Z≥0, k ∈ {0, . . . , Nl− 1}.
We write the system (22) compactly in the form ˙
x(t) = f (x(t)) + g(x(tlk)) + h(x(tl0)), (23) with
f (x) = col(DTTω, −M−1(DT∇U (ϕ) + Aω + Pd), 0, 0, 0),
g(x) = col(0, 0, Tb−1(Pg− Q−1b − Q−1c),
Tg−1(1λ − b − σ2ω + ρ11T(Pd− Pg)), τλ−11 T(P
d− Pg)),
t0 1 t02 t 0 k−1 t0k t 0 k+1 t11 t l k t0 N0= t 1 0 tlNl= t l+1 0 t0 0= 0 . . . .
power network dynamics (4) t
ISO-generator bidding process (Algorithm 1) ISO-generator bidding process power network dynamics frequency deviations power generation setpoints
Figure 1: Relation between time and iteration numbers in the time-triggered system (22). The lower time-axis corresponds to the continuous-time physical system (4) while the upper one corresponds to the time sequence {{tl
k} Nl
k=0} ∞
l=0 of the ISO-generator bidding process given in Algorithm 1. The arrows
pointing up indicate the frequency updates in the bidding dynamics while the arrows pointing down correspond to update of the power generation levels in the physical system. As indicated, for each l ∈ Z≥0the lower index k is reset once it reaches k = Nl∈ Z≥1, i.e., tlN
l= t
l+1
0 for all l ∈ Z≥0. With this notation, note that the continuous-time
dynam-ics (15) corresponds to ˙
x(t) = f (x(t)) + g(x(t)) + h(x(t)). (24) Since supϕ∈Rn−1k∇2U (ϕ)k < ∞ and g, h are linear, it follows that f, g, h are globally Lipschitz (we denote by Lf, Lg, Lh their Lipschitz constants, respectively). When
viewed as a continuous-time system, the dynamics (22) has a discontinuous right-hand side, and therefore we consider its solutions in the Carath´eodory sense, cf. [29].
B. Sufficient condition on triggering times for stability In this section we establish conditions on the time sequence that guarantee that the solutions of (22) are well-defined and retain the convergence properties of (15). Specifically, we determine a sufficient condition on the inter-sampling times tl
k+1− t l
k for bidding and t l+1 k − t
l
k for market clearing that
ensure local asymptotic convergence of (23) to the equilib-rium ¯x of the continuous-time system (15).
Our strategy to accomplish this relies on the robustness properties of (15) characterized in Theorem IV.1 and the fact that the time-triggered implementation, represented by (23), can be regarded as an approximation of the continuous-time dynamics, represented by (24). We use the Lyapunov function W defined by (18) and examine the mismatch between both
dynamics to derive upper bounds on the inter-event times that guarantee that W is strictly decreasing along the
time-triggered system (22).
Theorem V.1. (Local asymptotic stability of time-triggered implementation): Consider the time-triggered implementa-tion (22) of the interconnection between the ISO-generator bidding processes and the power network dynamics. With the notation of Theorem IV.1, let
¯ ξ := 1 Lf+ Lg log1 + β(Lf+ Lg) L(LWLh+ β) , (25) ¯ ζ := 1 Lf log1 + Lf(α − β) Lg(LLW + α) + (α − β)(Lf+ Lg) ,
where0 < β < α, L := Lf+Lg+Lh, andLW is the Lipschitz
constant of ∇W. Assume the time sequence {{tlk} Nl
k=0}∞l=0
satisfies, for someζ ∈ (0, ζ) and ξ ∈ (0, ξ),
ζ ≤ tl+10 − tl
0≤ ζ and ξ ≤ tlk− tlk−1≤ ξ, (26)
for alll ∈ Z≥0andk ∈ [Nl]. Then, ¯x is locally asymptotically
stable under(22).
We refer the reader to Appendix B for the proof of The-orem V.1. The uniform lower bounds ζ and ξ in (26) ensure that the solutions of the time-triggered implementation (22) are well-defined, avoiding Zeno behavior. Theorem V.1 implies that convergence is guaranteed for any constant stepsize imple-mentation, where the sufficiently small stepsize satisfies (26). However, the result of Theorem V.1 is more general and does not require constant stepsizes. Another interesting observation is that the upper bounds can be calculated without requiring any information about the equilibrium ¯x. This is desirable, as this equilibrium is not known beforehand and must be determined by the algorithm itself.
Remark V.2. (General cost functions, generator box con-straints, and elastic demand): We briefly discuss here the feasibility and challenges involved in extending the present treatment along different directions.
General cost functions: The results presented above also hold for general convex (instead of quadratic) cost func-tions. Specifically, if there exist c1, c2 ∈ R>0 such that
0 ≺ c1I ∇2C(Pg) c2I for all Pg ∈ Rn, then
one can establish LISS of the closed-loop system and compute bounds on the inter-event times to guarantee local asymptotic stability. However, this setting requires a significant more technical derivation and therefore it is omitted in the present work.
Generator box constraints: An interesting extension would be the inclusion in (7) of box constraints on power gen-eration. Such extension would require the re-examination of the existence and uniqueness of Nash equilibria, as generators may take advantage of them to arbitrarily increase their individual profit, and the modification of the proposed algorithm to accommodate these constraints, e.g., by using projections in the dynamics. The latter makes it difficult to asses LISS as the characterization of the input-to-state stability properties of projected dy-namics systems is still open.
Elastic demand: Our treatment here has considered inelastic loads. If instead loads are flexible, i.e., the demand is elas-tic, then consumers would also react to prices, possibly in a price-taking manner as in a Cournot-type of competitive market. This implies that the game-theoretic model has
to be revised to accommodate loads that change with generator bids/prices. The investigation of how to carry this out is an interesting topic for future research. •
VI. SIMULATIONS
In this section we illustrate the convergence properties of the interconnected time-triggered system (22). We consider the IEEE 14-bus power network depicted in Figure 2, where each node has one generator and one load according to model (2). We assume costs at each node i ∈ [14] of the form
Ci(Pgi) =12qiPgi2 + ciPgi,
with qi > 0 and ci ≥ 0. In the original IEEE 14-bus
bench-mark model, nodes 1, 2, 3, 6, 8 have synchronous generators while the other nodes are load nodes and have no power generation. We replicate this by suitably choosing the cost at the load nodes such that the optimizer of the economic dispatch problem (7) is zero at them. In addition, we choose Mi ∈ [4, 5.5] for generator nodes i ∈ {1, 2, 3, 6, 8} and
Mi 1 for the load nodes. We set Ai ∈ [1.5, 2.5], Vi ∈
[1, 1.06], Tbi ∈ [0.0005, 0.001], Tgi = 13.5 for all i ∈ [14]
and ρ = 900. The other parameter values for the ISO dynamics (12b) are τλ= 0.0004, ρ = 3, σ = 17. 1 2 3 4 5 6 7 8 9 10 11 12 13 14
Figure 2: Schematic of the modified IEEE 14-bus benchmark. Each edge represents a transmission line. Red nodes represent loads. All the other nodes represent synchronous generators, with different colors that match the ones used in Figures 3 and 5. The physical dynamics are modeled by (2).
At time t = 0 s, the inelastic load (in MW’s) is given by Pd= (0, 20, 86, 43, 7, 10, 0, 0, 27, 8, 3, 6, 12, 14).
Initially, we set (q1, q2, q3, q6, q8) = (22, 128, 45, 60, 30),
(c1, c2, c3, c6, c8) = (7.5, 7.5, 7.5, 7.5, 7.5) and qi =
1500, ci = 26 for the remaining nodes. The time-triggered
system (22) is initialized at steady state at the optimal gener-ation level
(Pg1, Pg2, Pg3, Pg6, Pg8) = (85, 15, 42, 31, 63),
and with Pgi = 0 for all other nodes. Figures 3-5 depict the
simulation of the time-triggered system for different triggering times. At t = 1 s all the loads are increased by 10% and we set ci = 28 for the load nodes. As observed in all figures,
the trajectories converge to a new efficient equilibrium with optimal power generation level
(Pg1, Pg2, Pg3, Pg6, Pg8) = (94, 16, 46, 34, 69),
and Pgi = 0 for all other nodes. Furthermore, at steady
state the generators all bid equal to the Lagrange multiplier which, by Proposition III.1, corresponds to an efficient Nash equilibrium.
At t = 15 s the cost functions of the generators are changed to (q1, q2, q3, q6, q8) = (23, 116, 48, 63, 38),
(c1, c2, c3, c6, c8) = (7.5, 6, 13.5, 15, 10.5) and qi =
1500, ci = 33 for the remaining nodes. As a result, the
optimal dispatch of power changes. Due to the changes of the power generation, a temporary frequency imbalance occurs. As illustrated in Figures 3-5, the power generations converge to the new optimal steady state given by
(Pg1, Pg2, Pg3, Pg6, Pg8) = (108, 23, 40, 28, 60).
In addition, we observe that after each change of either the load or the cost function, the frequency is stabilized and the bids converge to a new efficient Nash equilibrium. The fact that the frequency transients are better in Figures 3-4 (with inter-event times of maximal 2 ms for bidding and on average respectively 50 ms, 62.5 ms for market clearing) than in Figure 5 (with inter-event times of 2 ms for bidding and 160 ms for market clearing) is to be expected given the longer inter-event times in the second case. A slight increase in the inter-event times for Figure 5 in either bid updating or market clearing time result in an unstable system. Figure 6 illustrates the evolution of the interconnected system with the primary/secondary and tertiary control layers separated and its loss of efficiency compared to the proposed integrated design.
VII. CONCLUSIONS
This paper has studied the joint operation of the economic dispatch and frequency regulation layers, which are tradition-ally separated in the control of power networks. The starting point of our design was a continuous-time bid update scheme coupled with the frequency dynamics whose equilibrium cor-responds to an efficient Nash equilibrium and zero frequency deviation. Building on the identification of a novel LISS-Lyapunov function for this dynamics, we have characterized its robustness properties against additive disturbances. We have exploited the LISS-property to propose a provably correct multi-rate hybrid implementation that combines the iterative nature of the fast bid updates and the slower power setpoint updates with the continuous frequency network dynamics. Our results show that real-time iterative bidding can successfully be interconnected with frequency control to increase efficiency while retaining stability of the power system. Future work will incorporate elastic demand, generator bounds, and power flow constraints in the formulation. We also wish to explore more general bids including piecewise constant and linear supply functions. Finally, we aim to develop distributed and opportunistic self-triggered implementations of the proposed dynamics, and characterize the convergence properties of data-driven optimization algorithms.
APPENDIXA PROOF OFTHEOREMIV.1
We structure the proof of Theorem IV.1 in two separate parts, corresponding to the inequalities (20a) and (20b), re-spectively.
Time (s) 0 10 20 30 -0.015 -0.01 -0.005 0 0.005 0.01 Frequency deviation (Hz)
(a) Evolution of the frequency deviations. After each sudden supply-demand imbalance, frequency is restored to its nominal value.
Time (s) 0 10 20 30 0 20 40 60 80 100 120 Power generation (MW)
(b) Evolution of the nodal power generations. After each change in the network, the power generation quantities converge to the optimal values deter-mined by (7). Time (s) 0 10 20 30 26 28 30 32 34 36 Bid ($/MWh)
(c) Evolution of the bids and the Lagrange mul-tiplier (dashed black line). As shown, the bids converge to the unique efficient Nash equilibrium. Figure 3: Simulations of the interconnection between the iterative bidding mechanism and the power network dynamics modeled by the time-triggered system (22). The colors in the graph corresponds to the nodes as depicted in Figure 2. We choose identical inter-event times given by tl
k− tlk−1= 2 ms, tl0− t l−1 0 =
50 ms for all l ∈ Z≥1, k ∈ [25]. As expected, the time-triggered system is asymptotically stable for sufficiently fast updates.
Time (s) 0 10 20 30 -0.02 -0.01 0 0.01 Frequency deviation (Hz)
(a) Evolution of the frequency deviations.
Time (s) 0 10 20 30 0 20 40 60 80 100 120 Power generation (MW)
(b) Evolution of each power generation.
Time (s) 0 10 20 30 26 28 30 32 34 36 Bid ($/MWh)
(c) Evolution of the bids & Lagrange multiplier. Figure 4: Simulations of the time-triggered system (22) under time-varying step sizes. We choose the time between two consecutive bid iterations randomly between 0.5 ms ≤ tl
k− tlk−1≤ 2 ms, for all l ∈ Z≥1, k ∈ [Nl], and we choose the number of bid iterations Nl∈ Z before market clearing occurs randomly
in the interval [20, 80]. Since the step sizes are sufficiently small, and therefore the mismatch of the time-triggered system with its continuous-time variant, the performance is similar compared to Figure 3.
A. Positive definiteness of Lyapunov function W
Let ¯x be the equilibrium of (15) satisfying the hypothesis. We now prove the existence of constants c1, c2, 0 > 0 such
that (20a) holds, given the constants 1, 2, 3 > 0. The
Hessian of W(eq. (18)) is given by a block-diagonal matrix
∇2W
(x) = blockdiag(H1(ϕ), H2) with upper left block as
H1(ϕ) = ∇2U (ϕ) 01D†TM 01M DT†T M ,
and the lower right block is given by
H2= 1 σ2 Tb −02Tg 0 −02Tg Tg −03Tg1 0 −031TTg τλ .
We will now show that there exists sufficiently small 0 such
that H1(ϕ), H2 are both positive definite for all x ∈ Ω. To
this end, let us define the function
U(η) = D†
TDΓ cos(η)D TD†T
T (27)
and note that U(DTD†T
T ϕ) = ∇2U (ϕ), implying that 0 <
U(γ1) ≤ ∇2U (ϕ) ≤ ∇2U (0) =U(0) for all x ∈ Ω, see (19).
Consequently, forD := 01DT†M , we have
U(γ1) D DT M | {z } K1 ≤ H1(ϕ) ≤ U(0) D DT M | {z } K2 , ∀x ∈ Ω.
By considering the Schur complements, the matrices K1, H2
are shown to be positive definite by choosing 0 > 0
suffi-ciently small such that
U(γ1) − 2 021D † TM D †T T > 0, Tb− 20 2 2Tg> 0, τλ− 20231TTbTg(Tb− 2022Tg)−11 > 0. (28) Next we define
c1:= min{λmin(K1), λmin(H2)}, (29)
c2:= max{λmax(K2), λmax(H2)}, (30)
where λmin(A), λmax(A) denote the smallest and largest
eigenvalue of the matrix A ∈ Rn×n. Note that c1, c2 > 0
and the following holds
Time (s) 0 10 20 30 -0.015 -0.01 -0.005 0 0.005 0.01 Frequency deviation (Hz)
(a) Compared to Figure 3a, there are more os-cillations and a larger overshoot of the frequency deviations. Time (s) 0 10 20 30 0 20 40 60 80 100 120 Power generation (MW)
(b) Evolution of the power generations at each node. Time (s) 0 10 20 30 26 28 30 32 34 36 Bid ($/MWh)
(c) Evolution of the bids and the Lagrange multi-plier. Compared to Figure 3c, more oscillations in the bids occur.
Figure 5: Simulations of the time-triggered system (22). Here we consider the case tl
k− tlk−1= 2 ms, tl0− t l−1
0 = 160 ms for all l ∈ Z≥1, k ∈ [80]. The
scenario is the same as in Figure 3. In this case however, the interconnected time-triggered system is only marginally stable; a small increase in either of the inter-event times results in an unstable system.
Time (s) 0 10 20 30 -0.02 -0.01 0 0.01 Frequency deviation (Hz)
(a) Evolution of the frequency deviations. Com-pared to Figures 3a-4a, there are more oscillations in the frequency deviations.
Time (s) 0 10 20 30 0 20 40 60 80 100 120 Power generation (MW)
(b) Evolution of each power generation. After primary and secondary controllers are activated at t = 1 s, optimal power sharing is lost.
Time (s) 0 10 20 30 3500 4000 4500 5000 5500 6000 Total cost ($/h)
(c) Evolution of the total generation costs (in black) compared to the optimal values calculated by (7). Activation of primary/secondary control, and changes in the cost function result in a loss of efficiency.
Figure 6: Simulations of swing equations with the primary/secondary and tertiary control layers separated. At time t = 1 s, the load is increased as in Figure 3 and decentralized primary/secondary controllers are activated to regulate the frequency but, as a result, optimal power sharing is lost. At t = 14 s the tertiary control layer is activated by resetting the setpoints optimally. After the change of the cost functions at t = 15 s, economic optimality is temporary lost again until the next time the tertiary control layer is activated (typically in the order of minutes).
Note that since W(¯x) = 0, ∇W(¯x) = 0, we have
W(x) = W(x) − W(¯x) = (x − ¯x)TR01∇W((x − ¯x)τ + ¯x) − ∇W(¯x) dτ = (x − ¯x)TR1 0 R1 0 τ ∇ 2W ((x − ¯x)τ θ + ¯x)dτ dθ(x − ¯x).
Since Ω is convex, it follows that xτ θ + (1 − τ θ)¯x ∈ Ω for all τ, θ ∈ [0, 1], x ∈ Ω. Consequently, by (31) we have
c1I ≤ ∇2W(xτ θ + (1 − τ θ)¯x) ≤ c2I, ∀τ, θ ∈ [0, 1],
and ∀x ∈ Ω. SinceR01R01τ dθdτ = 1
2, inequality (20a) follows.
B. Dissipation inequality
Here we establish the inequality (20b). First we consider the case without disturbance, i.e., d = 0. Given the equilibrium ¯
x of (15), we define ˜x := x − ¯x and likewise ˜ϕ, ˜ω, ˜b, ˜Pg, ˜λ.
Then, the system (15) reads as ˙˜
ϕ = DTTω,˜
M ˙˜ω = −DT(∇U (ϕ) − ∇U ( ¯ϕ)) − A˜ω + ˜Pg,
Tb˙˜b = ˜Pg− Q−1˜b,
TgP˙˜g= 1˜λ − ˜b − ρ11TP˜g− σ2ω,˜
τλ˙˜λ = −1TP˜g.
In addition, note that W (eq. (18)) takes the form
W(x) = V (x) + V(x), (32) V(x) = 01ϕ˜TDT†M ω − 02 σ2 ˜b TT gP˜g−0σ23λ1˜ TT gP˜g. (33)
Next, we determine the time-derivative of the individual terms of the candidate Lyapunov function W.
(0): First, observe from (16) that ˙
V = −ωTAω −σ12(b − ¯b)
TQ−1(b − ¯b)
−σρ2(Pg− ¯Pg)T11T(Pg− ¯Pg).
(1): The time-derivative of the first term of V satisfies d dtϕ˜ TD† TM ω = ˜ω TM D†T T D T Tω˜
− ˜ϕTDT†DT(∇U (ϕ) − ∇U ( ¯ϕ)) − ˜ϕTDT†A˜ω + ˜ϕTDT†P˜g.
By exploiting D†TDT = I, the second term is rewritten as
where we used that ∇U (ϕ) − ∇U ( ¯ϕ) = U (ϕ)(ϕ − ¯ϕ) with U (ϕ) =R1
0 ∇
2U ((ϕ − ¯ϕ)θ + ¯ϕ)dθ. (34)
Since U (ϕ) ≥ U(1γ) = DT†DΓ cos(1γ)DTD†T
T (see eq.
(27)) for all x ∈ Ω, we obtain
d dtϕ˜ TD† TM ω ≤ ˜ω TM D†T T D T Tω − ˜˜ ϕ TU(1γ) ˜ϕT − ˜ϕTD†TA˜ω + ˜ϕTD†TP˜g.
(2): For the second term of V the following holds: d dt˜b TT gP˜g = ˜PgTTgbP˜g− ˜PgTTgbQ−1˜b + ˜bT1˜λ − ˜bT˜b − ρ˜bT 11TP˜g− σ2˜bTω,˜ where we define Tgb := TgTb−1.
(3): Similarly, by defining Tgλ:= TgTλ−1 we obtain d dtλ1˜ TT gP˜g= − ˜PgTTgλ11TP˜g+ n˜λ2− ˜λ1T˜b − ρn˜λ1TP˜ g− σ2λ1˜ Tω.˜
By combining the above calculations, we can show that the time-derivative of W satisfies
˙
W= ˙V + ˙V≤ 120(x − ¯x)TPTΞP(x − ¯x).
where Ξ is given by (35) (see page 11), P takes the form
P = 0 I 0 0 0 0 0 1 σI 0 0 0 0 0 1 σI 0 0 0 0 0 1 σ I 0 0 0 0 , andM := MD†TT DT T + DTDT†M,T := Tgλ11T + 11TTgλ.
Next, we will show that there exists 0, 1, 2, 3> 0 such that
Ξ is positive definite.
This can be done by successive use of the Schur com-plement. In particular, for A ∈ Rn×n
, B ∈ Rn×m, C ∈ Rm×m, β > 0, recall that βA B BT C > 0 ⇐⇒ C > 0 & βA − BC−1BT > 0.
For successively applying this result to Ξ, given by (35), let
us first fix 1, 3> 0. Then 2can be chosen sufficiently large
such that lower-right 3 × 3 block submatrix of Ξ is positive
definite. Then we can choose a 0 > 0 sufficiently small
such that (28) holds and Ξ > 0. Here, note that choosing
0 smaller does not affect the positive definiteness of the
lower-right 3 × 3 block submatrix of Ξ. By construction of
0, 1, 2, 3, there exist constants c1, c2∈ R>0such that (20a)
holds for all x ∈ Ω, see also Section A-A. In addition, for this choice of we have that Ξ> 0 and, as a result, there exists
ˆ
α := 120λmin(PTΞP) > 0 such that
(∇W(x))TF (x) ≤ − ˆαkx − ¯xk2
for all x ∈ Ω. Next, we consider the case when the disturbance is present. Let χ satisfy 0 < χ < ˆα/(LWkBk). Then, by
exploiting the Lipschitz property of ∇W,
(∇W(x))T(F (x) + Bd) ≤ − ˆαkx − ¯xk2+ ∇W(x))TBd
≤ − ˆαkx − ¯xk2+ LWkBkkx − ¯xkkdk
≤ −( ˆα − LWkBkχ)kx − ¯xk2= −αkx − ¯xk2
with α := ˆα − LWkBkχ > 0 and thus (20b) holds. This
concludes the proof of Theorem IV.1. APPENDIXB
PROOF OFTHEOREMV.1
Here we prove Theorem V.1. To do so, we rely on Gron-wall’s inequality, which in general allows to bound the evo-lution of continuous-time and discrete-time signals described by differential and difference equations, respectively. Given the hybrid nature of the time-triggered dynamics (22), we rely on a version of Gronwall’s inequality for hybrid systems developed in [30]. Adapted for our purposes, it states the following. Proposition B.1. (Generalized Gronwall’s inequality [30]): Let t 7→ y(t) ∈ R be a continuous signal, t 7→ p(t) ∈ R be a continuously differentiable signal, r := {rj}k−1j=0 be a
nonnegative sequence of real numbers,q ≥ 0 a constant, and E := {tj}k+1j=0, k ∈ Z≥0 be a sequence of times satisfying
tj < tj+1 for all j ∈ {0, . . . , k}. Suppose that for all
t ∈ [t0, tk+1], the elements y, p, and r satisfy
y(t) ≤ p(t) + qRt
t0y(s)ds +
Pi(t)−1
j=0 rjy(tj+1) (36)
withi(t) := max{i ∈ Z≥0: ti ≤ t, ti∈ E} for t < tk+1 and
i(tk+1) := k. Then,
y(t) ≤ p(t0)h(t0, t) +
Rt
t0h(s, t) ˙p(s)ds (37) for all t ∈ [t0, tk+1] where for all t0≤ s ≤ t ≤ tk+1,
h(s, t) := exp q(t − s) +Pi(t)−1
j=i(s)log(1 + rj). (38)
We are now ready to prove Theorem V.1. Proof of Theorem V.1: Let {{tl
k} Nl
k=0} ∞
l=0 be a sequence
of times satisfying the hypotheses. Consider a trajectory t 7→ x(t) of (22) with x(0) belonging to a neighborhood of ¯x. The definition of this neighborhood will show up later. Our proof strategy involves showing the monotonic decrease of the function W (cf. (18)) along this arbitrarily chosen trajectory.
Consider any t ∈ R≥0 such that t 6∈ {tlk} Nl
k=0for any l ∈ Z≥0
and x(t) ∈ Ω where Ω is defined by (19). With a slight abuse of notation let l and k ∈ {0, . . . , Nl− 1} be fixed such that
t ∈ (tlk, tlk+1). Then, using the expression of F (x) = f (x) + g(x) + h(x) given in (24), one can write the evolution of x at t for the considered trajectory as
˙
x(t) = F (x(t)) + g(x(tlk)) − g(x(t)) + h(x(tl0)) − h(x(t)).
(I) Dissipation inequality: Note that at t the evolution of W
is equal to the dot product between the gradient of W and
right-hand side of the above equation. Hence, we get ˙ W(x(t)) = ∇W(x(t))> F (x(t)) + g(x(tlk)) − g(x(t)) + h(x(tl0)) − h(x(t)). (39) From (20b), we have ∇W(x(t))>F (x(t)) ≤ −αkx(t) − ¯xk2.
Moreover, since maps ∇W, g, and h are globally Lipschitz
and ∇W(¯x) = 0, one has k∇W(x(t))k ≤ LWkx(t) − ¯xk,
kg(x(tl
k)) − g(x(t))k ≤ Lgkx(tlk) − x(t)k, and kh(x(t l 0)) −
Ξ= ω b/σ Pg/σ λ/σ ϕ ω 2 0A − 1M −2σI 0 −3σ1 1AD †T T b σ −2σI −22I + 2 0Q −1 − 2(Q−1Tgb+ ρ11T) (2− 3)1 0 Pg σ 0 −2(TgbQ −1 + ρ11T) 22Tgb+2 0ρ11 T − 3T −3nρ1 −1σDT†T λ σ −3σ1 T ( 2− 3)1T −3nρ1T 2n3 0 ϕ 1DT†A 0 −1σD†T 0 21U(1γ) (35)
h(x(t))k ≤ Lhkx(tl0) − x(t)k. Using these bounds in (39), we
get ˙ W(x(t)) ≤ −αkx(t) − ¯xk2+ LWkx(t) − ¯xk Lgkx(t) − x(tl k)k + Lhkx(t) − x(tl0)k . (40)
Next, we provide bounds on kx(t)−x(tl
k)k and kx(t)−x(tl0)k
in terms of kx(t) − ¯xk, t − tl
k, and t − tl0. To reduce the
notational burden, we drop the superscript l from the time instances {tli}Nl
i=1. In addition, we define
xk:= x(tk), ζk(t) := t − tk,
ζjk:= ζj(tk) = tk− tj, ξl(t) := ζ0(t) = t − t0.
(II) Bounds on kx(t) − x(tl
k)k: Note that x(t) can be written
using (23) as the line integral
x(t)−xk= Rt tkf (x(s))ds + ζk(t)g(xk) + ζk(t)h(x0) =Rt tk(f (x(s)) − f (xk))ds + ζk(t)(f (xk) − f (¯x)) + ζk(t)(g(xk) − g(¯x) + h(x0) − h(¯x)). (41)
Above, we have added and subtracted ζk(t)f (xk) and
sub-tracted f (¯x) + g(¯x) + h(¯x) as ¯x is an equilibrium. Using Lipschitz bounds and triangle inequality in (41) we obtain
kx(t) − xkk ≤ Lf
Rt
tkkx(s) − xkkds (42) + ζk(t)(Lf+ Lg)kxk− ¯xk + ζk(t)Lhkx0− ¯xk.
We wish to obtain an upper bound on kx(t) − xkk that is
independent of the state at times s ∈ (tk, t). To this end, we
employ Gronwall’s inequality as stated in its general form in Proposition B.1. Drawing a parallelism between the notations, for (42), we consider E = ∅, r = 0, y(t) = kx(t) − xkk, q =
Lf, p(t) = ζk(t)(Lf + Lg)kxk− ¯xk + ζk(t)Lhkx0− ¯xk and
replace t0 by tk. With these choices, the hypothesis (36) is
satisfied (as it exactly corresponds to (42) above). Then, with the notation of Proposition B.1, we have
h(s, t) = eLf(t−s), ˙ p(s) = (Lf + Lg)kxk− ¯xk + Lhkx0− ¯xk, p(tk) = 0, and (37) reads as kx(t) − xkk ≤ Rt tke Lf(t−s) (L f+Lg)kxk− ¯xk+Lhkx0− ¯xkds = 1 +Lg Lfkxk− ¯xk(e Lfζk(t)− 1) (43) +Lh Lfkx0− ¯xk(e Lfζk(t)− 1).
Bounding the above inequality using the triangle inequality collecting coefficients of kx(t) − xkk on the left-hand side,
and rearranging gives
kx(t) − xkk ≤ Lh(eLfζk(t)− 1) Lf− (Lf+ Lg)(eLfζk(t)− 1) kx0− ¯xk + (Lf+ Lg)(e Lfζk(t)− 1) Lf− (Lf+ Lg)(eLfζk(t)− 1) kx(t) − ¯xk. (44) (III) Bounds on kx(t) − x(tl
0)k: Our next step is to provide
an upper bound on the term kx(t) − x0k. Recall that the
considered trajectory satisfies (23) and so, the line integral over the interval [t0, t] gives
x(t) − x0= Rt t0f (x(s))ds + Pk−1 j=0ζ j+1 j g(xj) + ζk(t)g(xk) + ξl(t)h(x0).
As done before, on the right-hand side, we add and subtract the terms ξl(t)f (x0) and ξl(t)g(x0) and then subtract f (¯x) +
g(¯x) + h(¯x). This gives us x(t) − x0=R t t0(f (x(s)) − f (x0))ds +Pk−1 j=0ζ j+1 j (g(xj) − g(x0)) + ζk(t)(g(xk) − g(x0)) + ξl(t)(f (x0) − f (¯x) + g(x0) − g(¯x) + h(x0) − h(¯x)).
By defining L := Lf+ Lg+ Lh, taking the norms, using the
global Lipschitzness, we obtain from above
kx(t) − x0k ≤ Lf Rt t0kx(s) − x0kds + ξ l(t)Lkx 0− ¯xk + LgP k−1 j=0ζ j+1 j kxj− x0k + Lgζk(t)kxk− x0k.
Consider any ˆt ∈ [t, tk+1] and note that ζk(t) ≤ ζk(ˆt).
Using this bound and the fact that the first term in the above summation is zero, we write
kx(t)−x0k ≤ LfR t t0kx(s) − x0kds + ξ l(t)Lkx 0− ¯xk + LgPk−2j=0ζj+1j+2kxj+1− x0k + Lgζk(ˆt)kxk− x0k.
We now apply Proposition B.1 to give a bound for the left-hand side independent of x(s), s ∈ (t0, t]. In order to do
so, the elements corresponding to those in the Gronwall’s inequality are: E = {tj}kj=0∪ {ˆt}, y(t) = kx(t) − x0k, p(t) =
ξl(t)Lkx
0− ¯xk, q = Lf, rj = Lgζj+1j+2 for j = 0, . . . , k − 2,
and rk−1= ˆt − tk. From Proposition B.1, we get
kx(t) − x0k ≤ Lkx0− ¯xk Rt t0h(s, t)ds, (45) where h(s, t) = exp Rt sLfdT + Pk−2 j=i(s)log(1 + ζ j+2 j+1Lg) +
Using log(1 + x) ≤ x for x ≥ 0 and the fact that the exponential is a monotonically increasing function, we get
h(s, t) ≤ expLf(t − s) + LgPk−2j=i(s)ζj+1j+1+ Lgζk(ˆt)
.
By noting that s ≤ i(s) + 1 and t ≤ ˆt, we can upper bound the right-hand side as h(s, t) ≤ exp Lf(t − s) + Lg(ˆt − s).
Since ˆt was chosen arbitrarily in the interval [t, tk+1], we
pick it equal to t. Thus, h(s, t) ≤ exp ((Lg+ Lf)(t − s)).
Substituting this inequality in (45) yields kx(t) − x0k ≤ Lkx0− ¯xk Rt t0e (Lf+Lg)(t−s)ds = L L f+Lg(e (Lf+Lg)ξl(t)− 1)kx 0− ¯xk. (46)
This inequality when used in the right-hand side of the triangle inequality kx0− ¯xk ≤ kx(t) − x0k + kx(t) − ¯xk yields after
rearrangement the following
kx0− ¯xk ≤ Lf+ Lg Lf+ Lg− L(e(Lf+Lg)ξ l(t) − 1)kx(t) − ¯xk. (47) Subsequently, using the above bound in (46) gives
kx(t) − x0k ≤
L(e(Lf+Lg)ξl(t)− 1) Lf+ Lg− L(e(Lf+Lg)ξl(t)− 1)
kx(t) − ¯xk.
(48) Combining inequalities (44) and (47) we obtain
kx(t) − xkk ≤ Lh(eLfζk(t)− 1) Lf− (Lf+ Lg)(eLfζk(t)− 1) · Lf+ Lg Lf+ Lg− L(e(Lf+Lg)ξ l(t) − 1)kx(t) − ¯xk + (Lf+ Lg)(e Lfζk(t)− 1) Lf− (Lf+ Lg)(eLfζk(t)− 1) kx(t) − ¯xk. (49)
(IV) Monotonic decrease of W: Note first that
follow-ing (48) and usfollow-ing the bound ξl(t) ≤ ¯ξ yields
kx(t) − x0k ≤
L(e(Lf+Lg) ¯ξ− 1) Lf+ Lg− L(e(Lf+Lg) ¯ξ− 1)
kx(t) − ¯xk.
Using the definition of ¯ξ, one gets e(Lf+Lg) ¯ξ − 1 =
β(Lf+Lg)
L(LWLh+β). Substituting this value in the above inequality and simplifying the expression provides us
kx(t) − x0k ≤ (β/(LWLh))kx(t) − ¯xk. (50)
In a similar way, using the bound on ξl(t) and substituting the
value of e(Lf+Lg) ¯ξ− 1 in (49) and simplifying yields
kx(t) − xkk ≤
(eLfζk(t)− 1)(L + β/L
W)
Lf− (Lf+ Lg)(eLfζk(t)− 1)
kx(t) − ¯xk.
Note that ζk(t) ≤ ¯ζ. Using this bound and the definition of ¯ζ
in the above inequality gives kx(t) − xkk ≤ Lα−β
WLgkx(t) − ¯xk. (51) Finally, substituting (50) and (51) in (40) and using the fact that β < α, we obtain ˙W(x(t)) < 0. Recall that t ∈ R≥0
was chosen arbitrarily satisfying t 6∈ {tlk}Nl
k=1for any l ∈ Z≥0.
Therefore, W monotonically decreases at all times along the
trajectory except for a countable number of points. Further, the map t 7→ W(x(t)) is continuous. Therefore, we conclude that
the trajectory initialized in a compact level set of Wcontained
in Ω converges asymptotically to the equilibrium point ¯x. This
completes the proof.
REFERENCES
[1] T. W. Stegink, A. Cherukuri, C. De Persis, A. J. van der Schaft, and J. Cort´es, “Integrating iterative bidding in electricity markets and frequency regulation,” in American Control Conference, Milwaukee, Wisconsin, USA, 2018, to appear.
[2] B.-M. Hodge and M. Milligan, “Wind power forecasting error distribu-tions over multiple timescales,” in Power and Energy Society General Meeting. IEEE, 2011, pp. 1–8.
[3] S.-H. Huang, D. Maggio, K. McIntyre, V. Betanabhatla, J. Dumas, and J. Adams, “Impact of wind generation on system operations in the deregulated environment: Ercot experience,” in Power & Energy Society General Meeting. IEEE, 2009, pp. 1–8.
[4] F. D¨orfler, J. W. Simpson-Porco, and F. Bullo, “Breaking the hierar-chy: Distributed control and economic optimality in microgrids,” IEEE Transactions on Control of Network Systems, vol. 3, no. 3, pp. 241–253, 2016.
[5] X. Zhang and A. Papachristodoulou, “A real-time control framework for smart power networks: Design methodology and stability,” Automatica, vol. 58, pp. 43–50, 2015.
[6] S. Trip, M. B¨urger, and C. De Persis, “An internal model approach to (optimal) frequency regulation in power grids with time-varying voltages,” Automatica, vol. 64, pp. 240–253, 2016.
[7] N. Li, C. Zhao, and L. Chen, “Connecting automatic generation control and economic dispatch from an optimization view,” IEEE Transactions on Control of Network Systems, vol. 3, no. 3, pp. 254–264, 2016. [8] Y. Zhang, M. Hong, E. Dall’Anese, S. Dhople, and Z. Xu, “Distributed
controllers seeking AC optimal power flow solutions using ADMM,” IEEE Transactions on Smart Grid, 2018, to appear.
[9] R. Patel, C. Li, X. Yu, and B. McGrath, “Optimal automatic generation control of an interconnected power system under network constraints,” IEEE Transactions on Industrial Electronics, vol. 65, no. 9, pp. 7220– 7228, 2018.
[10] C. Zhao, E. Mallada, S. Low, and J. Bialek, “A unified framework for frequency control and congestion management,” in Power Systems Computation Conference (PSCC), 2016. IEEE, 2016, pp. 1–7. [11] Z. Wang, F. Liu, S. H. Low, C. Zhao, and S. Mei, “Distributed frequency
control with operational constraints, part II: Network power balance,” IEEE Transactions on Smart Grid, pp. 1–1, 2017.
[12] F. L. Alvarado, J. Meng, C. L. DeMarco, and W. S. Mota, “Stability analysis of interconnected power systems coupled with market dynam-ics,” IEEE Transactions on Power Systems, vol. 16, no. 4, pp. 695–701, 2001.
[13] D. J. Shiltz, M. Cvetkovi´c, and A. M. Annaswamy, “An integrated dynamic market mechanism for real-time markets and frequency reg-ulation,” IEEE Transactions on Sustainable Energy, vol. 7, no. 2, pp. 875–885, 2016.
[14] D. J. Shiltz, S. Baros, M. Cvetkovi´c, and A. M. Annaswamy, “Integration of automatic generation control and demand response via a dynamic regulation market mechanism,” IEEE Transactions on Control Systems Technology, 2018, to appear.
[15] T. W. Stegink, C. De Persis, and A. J. van der Schaft, “A unifying energy-based approach to stability of power grids with market dynamics,” IEEE Transactions on Automatic Control, vol. 62, no. 6, pp. 2612–2622, 2017. [16] R. Singh, P. R. Kumar, and L. Xie, “The ISO problem: Decentralized stochastic control via bidding schemes,” in 53rd Annual Allerton Con-ference on Communication, Control, and Computing (Allerton),. IEEE, 2015, pp. 552–559.
[17] R. Singh, P. R. Kumar, and L. Xie, “Decentralized control via dynamic stochastic prices: The independent system operator problem,” IEEE Transactions on Automatic Control, 2018, publication date: 30 January. [18] A. Cherukuri and J. Cort´es, “Iterative bidding in electricity markets: rationality and robustness,” arXiv preprint arXiv:1702.06505, 2017, submitted to IEEE Transactions on Control of Network Systems. [19] T. W. Stegink, A. Cherukuri, C. De Persis, A. J. van der Schaft, and
J. Cort´es, “Frequency-driven market mechanisms for optimal dispatch in power networks,” IEEE Transactions on Automatic Control, 2018, submitted.
[20] J. Machowski, J. W. Bialek, and J. R. Bumby, Power System Dynamics: Stability and Control, 2nd ed. Ltd: John Wiley & Sons, 2008. [21] P. Kundur, Power System Stability and Control. Mc-Graw-Hill
Engi-neering, 1993.
[22] F. D¨orfler and F. Bullo, “Synchronization in complex networks of phase oscillators: A survey,” Automatica, vol. 50, no. 6, pp. 1539–1564, 2014. [23] A. Cherukuri and J. Cort´es, “Decentralized Nash equilibrium learning by strategic generators for economic dispatch,” in American Control Conference. IEEE, 2016, pp. 1082–1087.
[24] D. Fudenberg and J. Tirole, Game Theory. Cambridge, MA: MIT Press, 1991.
[25] S. Stoft, “Power system economics,” Journal of Energy Literature, vol. 8, pp. 94–99, 2002.
[26] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011.
[27] E. D. Sontag and Y. Wang, “New characterizations of input-to-state stability,” IEEE Transactions on Automatic Control, vol. 41, no. 9, pp. 1283–1294, 1996.
[28] H. K. Khalil, Nonlinear systems. Prentice Hall, New Jersey, 1996, vol. 3.
[29] J. Cort´es, “Discontinuous dynamical systems - a tutorial on solutions, nonsmooth analysis, and stability,” IEEE Control Systems Magazine, vol. 28, no. 3, pp. 36–73, 2008.
[30] N. Noroozi, D. Neˇsi´c, and A. R. Teel, “Gronwall inequality for hybrid systems,” Automatica, vol. 50, no. 10, pp. 2718–2722, 2014.
Tjerk Stegink is a Ph.D. candidate at the Engineer-ing and Technology Institute, Faculty of Science and Engineering, University of Groningen, the Nether-lands. He received his B.Sc. (2012) in Applied Math-ematics and M.Sc. (2014, cum laude) in Applied Mathematics (Systems, Control and Optimization) from the same university. His main research interests are in distributed optimization and control of power systems.
Ashish Cherukuri is a Postdoctoral Researcher at the Automatic Control Laboratory, ETH Z¨urich, Switzerland. He received the B.Tech degree in me-chanical engineering from the Indian Institute of Technology, Delhi, India, in 2008, the M.Sc. de-gree in mechanical engineering from ETH, Zurich, Switzerland, in 2010, and the Ph.D. degree in me-chanical and aerospace engineering at the Univer-sity of California, San Diego, CA, USA, in 2017. His research interests include dynamical systems, game theory, distributed optimization, and power networks. He received the Excellence Scholarship in 2008 from ETH Z¨urich, the Focht-Powell Fellowship in 2012 from the University of California, San Diego, and the Robert E. Skelton Systems and Control Dissertation Award in 2017 from the University of California, San Diego. He was a finalist for Best Student Paper Award in the American Control Conference, 2015.
Claudio De Persis is a Professor at the Engineering and Technology Institute, Faculty of Science and Engineering, University of Groningen, the Nether-lands. He received the Laurea degree in Electronic Engineering in 1996 and the Ph.D. degree in Sys-tem Engineering in 2000 both from the University of Rome “La Sapienza”, Italy. Previously he held faculty positions at the Department of Mechanical Automation and Mechatronics, University of Twente and the Department of Computer, Control, and Management Engineering, University of Rome “La Sapienza”. He was a Research Associate at the Department of Systems Science and Mathematics, Washington University, St. Louis, MO, USA, in 2000-2001, and at the Department of Electrical Engineering, Yale University, New Haven, CT, USA, in 2001-2002. His research interest is focused on automatic control. He was an Editor of the International Journal of Robust and Nonlinear Control (2006-2013), an Associate Editor of the IEEE Transactions on Control Systems Technology (2010-2015), and of the IEEE Transactions on Automatic Control (2012-2015). He is currently an Associate Editor of Automatica (2013-present) and of IEEE Control Systems Letters (2017-present).
Arjan van der Schaft received the undergraduate and PhD degrees in mathematics from the University of Groningen, the Netherlands. In 1982 he joined the Department of Applied Mathematics, University of Twente, where he was appointed as full professor in Mathematical Systems and Control Theory in 2000. In September 2005 he returned to his Alma Mater as professor in Mathematics. Arjan van der Schaft is Fellow of the Institute of Electrical and Electronics Engineers (IEEE), Fellow of the International Fed-eration of Automatic Control (IFAC), and was the 2013 recipient of the 3-yearly awarded Certificate of Excellent Achievements of the IFAC Technical Committee on Nonlinear Systems. He was Invited Speaker at the International Congress of Mathematicians, Madrid, 2006. Books authored by him include Variational and Hamiltonian Control Systems (1987, with P.E. Crouch), Nonlinear Dynamical Control Systems (1990, 2016, with H. Nijmeijer), L2-Gain and Passivity Techniques in Nonlinear Control (1996, 2000, 2017), An Introduction to Hybrid Dynamical Systems (2000, with J.M. Schumacher), and Port-Hamiltonian Systems: An Introductory Overview (2014, with D. Jeltsema).
Jorge Cort´es (M’02-SM’06-F’14) received the Li-cenciatura degree in mathematics from Universidad de Zaragoza, Zaragoza, Spain, in 1997, and the Ph.D. degree in engineering mathematics from Uni-versidad Carlos III de Madrid, Madrid, Spain, in 2001. He held postdoctoral positions with the Uni-versity of Twente, Twente, The Netherlands, and the University of Illinois at Urbana-Champaign, Urbana, IL, USA. He was an Assistant Professor with the Department of Applied Mathematics and Statistics, University of California, Santa Cruz, CA, USA, from 2004 to 2007. He is currently a Professor in the Department of Mechanical and Aerospace Engineering, University of California, San Diego, CA, USA. He is the author of Geometric, Control and Numerical Aspects of Nonholonomic Systems (Springer-Verlag, 2002) and co-author (together with F. Bullo and S. Mart´ınez) of Distributed Control of Robotic Networks (Princeton University Press, 2009). He has been an IEEE Control Systems Society Distinguished Lecturer (2010-2014) and is an elected member for 2018-2020 of the Board of Governors of the IEEE Control Systems Society. His current research interests include distributed control and optimization, network science, opportunistic state-triggered control and coordination, reasoning under uncertainty, and distributed decision making in power networks, robotics, and transportation.
