Optimal secondary frequency regulation with on-off loads in power networks

Load side participation can provide support to the power network by appropriately adapting the demand when required. In addition, it enables an economically improved power allocation. In this study, we consider the problem of providing an optimal power allocation among generation and on-off loads within the secondary frequency control timeframe. In particular, we consider a mixed integer optimization problem which ensures that the secondary frequency control objectives (i.e. generation/demand balance and the frequency attaining its nominal value at steady state) are satisfied. We present analytical conditions on the generation and on-off load profiles such that an $\epsilon$-optimality interpretation of the steady state power allocation is obtained, providing a non-conservative value for $\epsilon$. Moreover, we develop a hierarchical control scheme that provides on-off load values that satisfy the proposed conditions. We study the interaction of the proposed control scheme with the physical dynamics of the power network and provide analytic stability guarantees. Our results are verified with numerical simulations on the Northeast Power Coordinating Council (NPCC) 140-bus system, where it is demonstrated that the proposed algorithm yields a close to optimal power allocation.


I. INTRODUCTION
Motivation and literature review: The penetration of renewable sources of generation in power networks is expected to grow over the next years, driven by technological advances and environmental concerns [2], [3]. This will make generation more intermittent, resulting in more frequent generationdemand imbalances that may harm power quality and even cause blackouts. In addition, a large penetration of renewable generation will limit the controllability of generation [4]. Hence, novel challenges need to be overcome to enable the safe and efficient operation of power networks, motivating the analytical study of their stability properties.
Demand side participation is considered by many to be a key way to address the above problem [5], due to the ability of loads to provide a fast response when required. In addition, load side participation offers advantages such as lower fuel consumption and greenhouse gas emissions, and better localization of disturbances. However, the introduction A preliminary version of this work has appeared in [1]. This manuscript provides additional results and the analytic proofs of the main results. Moreover, it includes additional discussion and simulations that demonstrate the impact of the proposed analysis. of a large number of active elements in the power grid raises the issue of fairness in power allocation, which can be interpreted as a problem of economic optimality. In addition, it motivates a transition from, traditionally implemented, centralized control schemes towards hierarchical, distributed and decentralized schemes. Such schemes allow a limited share of information among agents, which can improve cybersecurity and reduce the expense of the necessary communication infrastructure. In addition, they have the potential to respect the privacy of sensitive information, such as measurements and cost functions. Such schemes have the potential to yield scalable designs, where controller parameters do not need to be recalibrated when generation or controllable demand is introduced or removed from the power grid, enabling plug and play capabilities.
The above have motivated several studies to consider distributed and decentralized schemes for generation and demand control in power networks within the primary and secondary frequency control timeframes, where the objectives are to ensure generation-demand balance and that frequency attains its nominal value at steady state respectively. In particular, [6] proposed a decentralized scheme for demand control within the primary frequency control timeframe. In addition, [7] considered the stability and economic optimality of controllable demand schemes. These ideas were extended in [8], [9], which considered the stability and optimality properties of wide classes of generation and controllable demand schemes. Moreover, [10] proposed decentralized and distributed conditions that guaranteed the stability of a broad class of dynamics for generation and demand.
The stability and optimality properties of distributed schemes within the secondary frequency control timeframe were studied in [11], [12], [13], [14], [15]. These properties were demonstrated for a rich class of generation and demand dynamics in [16], [17], by considering passivity-related properties. In addition, [18], [19], [20] considered distributed schemes that also take into account various static or transient operational constraints of the power network. Moreover, [21] considered the effect of time-delays on the stability properties of a distributed optimality scheme. To achieve optimality, these studies construct appropriate optimization problems that enable an economically optimal power allocation and ensure that system equilibria coincide with the solutions to these problems. Furthermore, [22] proposed a decentralized scheme for optimal power sharing by considering a linear quadratic regulator problem. A thorough survey on distributed optimization schemes in power networks can be found in [23].
On many occasions, loads are described by on and off states and hence a continuous representation cannot accurately characterize their behavior. The possible switching of loads has been studied in [24], [25], which considered the problem of using on-off loads for ancillary support in power networks and provided stability guarantees for arbitrary network topologies within the primary and secondary frequency control timeframes respectively. In addition, [24] enabled an optimality interpretation of the resulting equilibria by means of a centralized information structure. Loads that switch into different modes of operation, corresponding to nominal and urgent circumstances respectively, were considered in [26]. Several studies considered the temperature-dependent, on-off behavior of loads [27], [28], [29], proposing various control schemes for improved performance. In summary, considering on-off load behavior and improving the stability and optimality properties of power networks when those are incorporated is of high importance.
Contribution: This study considers the optimality and stability properties of the power network when controllable on-off loads are incorporated within the secondary frequency control timeframe. The discontinuous nature of on-off loads introduces several challenges, requiring tools from switching system analysis and mixed integer programming.
In particular, we formulate a mixed integer optimization problem that associates the cost of generation and on-off controllable loads such that the secondary frequency control objectives are attained, i.e. generation-demand balance is achieved and the frequency takes its nominal value at steady state. The binary nature of on-off loads makes the problem combinatorial and hence the task of finding an exact solution prohibitively difficult, particularly when the number of loads is large. We hence adopt an ǫ-optimality approach, where we aim to achieve a cost of no more than ǫ from the global minimum. We propose equilibrium conditions that, when satisfied, ensure that the steady state power allocation is ǫ-optimal and provide a non-conservative value for ǫ. Moreover, we propose a hierarchical control policy that yields an equilibrium allocation that satisfies these conditions, by determining the on-off load values. A significant feature of the proposed scheme is its ability to accommodate the addition and removal of onoff loads without recalibration. Furthermore, we consider a suitably adjusted version of a distributed 'Primal-Dual' scheme that has been widely used in the literature, in studies aiming to obtain an optimal power allocation among continuous generation and demand units. The interaction of these schemes with the physical dynamics of the power network renders their combination a switching system. For the combined system, we explain that no chattering behavior occurs and provide asymptotic stability guarantees. In particular, our results ensure the convergence of solutions to an equilibrium that satisfies the secondary frequency control objectives and is ǫ-optimal to the considered optimization problem. The proposed schemes and corresponding analysis are applicable to arbitrary (connected) network configurations.
Our analytic results are verified with numerical simulations on the NPCC 140-bus system, that demonstrate the stability and optimality properties of the proposed algorithm on a realistic setting. In particular, the simulation results demonstrate that the proposed scheme enables an ǫ-optimality interpretation of the steady state power allocation.
To the authors best knowledge this is the first study that: (i) Offers analytic equilibrium conditions that ensure an ǫoptimal power allocation among generation and on-off controllable loads within the secondary frequency control timeframe, providing a non-conservative value for ǫ. (ii) Proposes a scalable hierarchical control scheme that determines on-off load values such that stability is guaranteed and secondary frequency control objectives and an ǫ-optimal power allocation are attained at steady state.
Paper structure: In Section II we present the considered model for the power network and the problem statement. In Section III we present our optimality analysis and the proposed hierarchical control scheme to solve the considered optimization problem. In Section IV we study the interaction of the proposed scheme with the physical network and present our main stability and optimality results. Our analytical results are demonstrated with numerical simulations in Section V and conclusions are drawn in Section VI. Proofs of the main results are provided in the Appendix.
Notation: Real and integer numbers are denoted by R and Z respectively. The set of n-dimensional vectors with real and integer entries are denoted by R n and Z n respectively. The sets of positive and non-negative real numbers are denoted by R >0 and R ≥0 respectively. We use 0 n and 1 n to denote n-dimensional vectors with all elements equal to 0 and 1 respectively. The image of a vector x is denoted by Im(x). The cardinality of a discrete set S is denoted by |S|. The convex closure of a set A is denoted byĀ. Finally, a function f : R → R is said to be monotonically increasing (respectively decreasing) if for all x and y such that x ≤ y it holds that f (x) ≤ f (y) (respectively f (x) ≥ f (y)).

A. Network model
We describe the power network by a connected graph (N , E) where N = {1, 2, . . . , |N |} is the set of buses and E ⊆ N ×N the set of transmission lines connecting the buses. Furthermore, we use (i, j) to denote the link connecting buses i and j and assume that the graph (N , E) is directed with an arbitrary orientation, so that if (i, j) ∈ E then (j, i) / ∈ E. It should be noted that the form of the presented dynamics is unaltered by any change in the graph ordering, and all of our results are independent of the choice of direction. For each j ∈ N , we use N p j = {k : (k, j) ∈ E} and N s j = {k : (j, k) ∈ E} to denote the sets of buses that are predecessors and successors of bus j respectively. The following assumptions are made about the network: 1) Bus voltage magnitudes are |V j | = 1 p.u. for all j ∈ N . 2) Lines (i, j) ∈ E are lossless and characterized by the magnitudes of their susceptances B ij = B ji > 0.
3) Reactive power flows do not affect bus voltage phase angles and frequencies. 4) Relative phase angles are sufficiently small such that the approximation sin η ij = η ij is valid.
The first three conditions are widely used in the literature [7], [8], [16], in studies associated with frequency regulation. These assumptions are valid in medium to high voltage transmission systems since transmission lines are dominantly inductive and voltage variations are small. The fourth condition is justified from the fact that the relative phase angles among buses are small in nominal operating conditions. Note that the theoretical results presented in this paper are validated with numerical simulations in Section V, which consider a detailed model of the power system with realistic parameter values.
We use the swing equations to describe the rate of change of frequency at each bus (e.g. [30], [31]). In particular, at each bus we consider some mechanical power injection from generation, a set of controllable on-off loads, and uncontrollable frequency-dependent and frequency-independent demand. This motivates the following system dynamics, (1c) In system (1), variables p M j and ω j represent the mechanical power injection and the deviation from the nominal value 1 of the frequency at bus j respectively. Variable d c l,j represents the demand of controllable load l at bus j. The set of controllable loads at bus j is denoted by L j . Furthermore, we define the setL := {(l, j) : l ∈ L j , j ∈ N }, such that all pairs l ∈ L j , j ∈ N satisfy (l, j) ∈L. The variable d u j represents the uncontrollable frequency-dependent load and generation damping present at bus j. Furthermore, variables η ij and p ij represent the power angle difference, and the power transmitted from bus i to bus j respectively. The constant M j > 0 denotes the generator inertia. Moreover, the constant p L j denotes the frequency-independent load at bus j, and ℓ = 1 T |N | p L its aggregate value throughout the network. Within the rest of the manuscript, we let x * denote the equilibrium value of state x.

B. Generation and demand dynamics
We consider generation and frequency-dependent uncontrollable demand and frequency damping dynamics described by where γ j > 0, j ∈ N , are time constants, p c j is a local power command variable available for design (see Section III-A), and A j > 0 and κ j > 0, j ∈ N , are damping and droop coefficients respectively. Note that the analysis carried in this paper is valid for more general generation and demand dynamics, including cases of nonlinear and higher order dynamics, provided certain input-output conditions hold, following the analysis in e.g. [8], [10], [16], [25]. In this paper, we consider first-order generation and static uncontrollable demand dynamics for simplicity and to retain the focus of the paper on on-off loads.
We consider on-off loads described by where d l,j ∈ R + denotes the magnitude of load (l, j) ∈L.
Variable σ l,j ∈ B = {0, 1} denotes the switching state of the lth load at bus j. The dynamics of σ are discussed in Sections III and IV. Furthermore, we define the constants ρ l,j ∈ B, which denote the desired switching state for each load (l, j) ∈ L, selected by its user.

C. Optimal generation and on-off load control
In this subsection we consider how generation and on-off controllable loads should be adjusted such that the cost of their joint operation is minimized while at the same time the generation and demand are balanced. In particular, we let 1 2 q j (p M j ) 2 be the cost incurred when the generation is p M j . Furthermore, we let a cost c l,j be incurred when the switching state σ l,j is different than the desired state ρ l,j at some on-off load (l, j) ∈L. The cost function for on-off loads is given by We then consider the following optimization problem, called the hybrid optimal supply control problem (H-OSC), H-OSC: The equality constraint in (4) requires all the frequencyindependent demand to be matched by the total generation and on-off controllable demand. This ensures that when system (1)-(2) is at equilibrium, the frequency will be at its nominal value. The latter follows by summing (1b) at steady state over all j ∈ N and noting the equality constraint in (4), which using (1a) and (2b) at steady state results to The second constraint reflects that controllable loads take discrete values. The latter makes (4) a combinatorial problem and hence difficult to solve when a large number of on-off loads is considered.

D. Problem statement
Below we state the main problem we aim to solve. Problem 1: Design a control scheme for generation and onoff loads, described by (2a) and (3), that: (i) Allows on-off loads to be added or removed from the network without requiring recalibration. (ii) Enables decentralized stability guarantees.
(iii) Is applicable to any (connected) network topology. (iv) Ensures that the frequency attains its nominal value at steady state. (v) Provides an optimality interpretation of the steady state power allocation. The first condition requires a control scheme that does not need to be modified when on-off loads are added or removed from the network, something that is expected to frequently occur due to their large numbers. The second requires that the control scheme enables locally verifiable stability guarantees. The third condition requires that the control scheme is applicable to any connected network configuration, i.e. its parameters do not depend on the topology of the power network. Condition (iv) is the main objective of secondary frequency control, i.e. to ensure that the frequency takes its nominal value at equilibrium. The last condition aims to ensure that the incurred cost is close to the global minimum of the H-OSC problem (4).

A. Power Command Dynamics
In this section we consider the design of the dynamics for power command variables, used as inputs to (2a), which contribute to enable an optimal power allocation. We adopt a suitably adapted version of a scheme that has been widely used in the literature for distributed optimal secondary frequency regulation [16], [18], [19], usually referred as the 'Primal-Dual' scheme. In particular, we consider a communication network described by a connected graph (N ,Ẽ), whereẼ represents the set of communication lines among the buses. In addition, for each j ∈ N , we useÑ p j = {k : (k, j) ∈Ẽ} andÑ s j = {k : (j, k) ∈Ẽ} to denote the sets of buses that are predecessors and successors of bus j within the communication network respectively. We consider the following dynamics for the power command signal p c j , where τ j and τ ij are positive constants, p c i and p c j are the power command signals which are shared between communicating buses i and j, while variable ψ ij is a state of the controller that integrates the power command difference between communicating buses i and j. We note that the set of communication linesẼ and power lines E can be the same or different.
The dynamics in (5) achieve both the consensus of the communicated variable p c , something that can be exploited to guarantee optimality of the equilibrium point reached, and also that frequency attains its nominal value at steady state. These are analytically shown in Lemma 1 below, proven in the Appendix, which follows by considering (1a), (2b), (5a) and the summation of (1b) and (5b) over all j ∈ N at steady state. Note that since the dynamics for σ have not been defined yet, the lemma is stated for the case where σ takes any constant value in B |L| . Lemma 1: Let σ take any constant value in B |L| . Then, all equilibria of (1), (2), (3), (5) satisfy ω * = 0 |N | and p c, * ∈ Im(1 |N | ).
Remark 1: An issue of implementability may be raised due to the requirement for demand measurements in the Primal-Dual scheme (5). This issue have been extensively considered in the literature [16], [32], [19], [33], [34], where schemes with different information structures have been proposed to relax the load measurements requirement. Note that exact knowledge of the demand is not required for the stability properties of the power system, presented below, to hold. In addition, demand estimates may be obtained from historical data and frequency measurements. However, an inaccurate estimate of the demand will affect the equilibrium properties of the power system and might result in a suboptimal power allocation. Hence, there exists a trade-off between the accuracy in knowledge of the demand and optimality. Nevertheless, a reasonable demand estimate can be easily obtained, which will allow a close to optimal allocation.

B. Optimality analysis
The H-OSC problem (4) aims to optimize the power allocation between generation and on-off loads. The considered problem is combinatorial and hence difficult to analytically solve when a large number of on-off loads is considered. Hence, an aim of this paper is to enable a steady state power allocation cost that is close to the global minimum of the H-OSC problem (4). Below, we define the notion of an ǫ-optimal point that is used throughout the rest of this manuscript.
Definition 1: Given a cost function C f : Below, we provide conditions that, when satisfied, enable an optimality interpretation of the steady state power allocation among generation and on-off loads. For convenience, we let β = max (l,j)∈L d l,j be the largest on-off load and K = j∈N k j the aggregate droop respectively within the power network.
Proposition 1: Consider an equilibrium of (1), (2), (3), (5). If k j = q −1 j , j ∈ N holds and there exists ζ ∈ R such that of (4). The parameter ζ associates the equilibrium values of power command variables, p c, * , and the load costs per unit demand c l,j /d l,j . The difference between ζ and p c, * is responsible for the additional cost, which accounts for the ǫoptimality interpretation. In addition, the parameter θ defines the margin of possible difference between ζ and p c, * and the value of ǫ. This suggests that when ζ = p c, * , i.e. when θ = 0, then a globally optimum solution can be obtained. However, there is no guarantee that such an equilibrium is always feasible. The following lemma, proven in the Appendix, demonstrates that when θ = β, then an equilibrium that satisfies (6) always exists and hence ǫ-optimality can be guaranteed.

C. Hierarchical control scheme for on-off loads
In this section we present a hierarchical control scheme which aims to obtain a vector σ which enables an ǫ-optimality interpretation of the resulting equilibria. The scheme is described by Algorithm 1 and a schematic representation of its information flow is depicted in Figure 1. Below, we provide additional explanations and intuition on its implementation.
Algorithm parameters: We define γ l,j = c l,j /d l,j , (l, j) ∈ L, which denotes the cost per unit demand for load (l, j), and the set G = ∪ (l,j)∈L {−γ l,j , γ l,j }. In addition, we letβ = β+δ, where the positive constant δ is allowed to be arbitrarily small and must satisfy δ ∈ (0, min(min γi,j ,γ k,l ∈G |γ i,j − γ k,l |), min (l,j)∈L d l,j ). The parameterβ is closely associated with Proposition 1 since the proposed algorithm aims to satisfy (6) with θ =β. It should be noted that all references to (6) within this section imply its satisfaction with θ =β. In addition, we defineβ = β + δ/2, which is used within the proposed algorithm.
In addition, we define perturbed costsc = c + r, where the vector r ∈ (0, δ/2) |L| is assumed to contain no duplicate elements. The vectorc serves to resolve complicacies resulting when loads have identical costs per unit demand c l,j /d l,j , by enabling distinct perturbed costs.
Inputs, states, output and information flow: The proposed scheme assumes knowledge of parameter K at the coordinating controller, p L j at each bus j ∈ N and d l,j ,c l,j at each load Algorithm 1: Hierarchical optimality scheme.
(l, j) ∈L, which can be regarded as its inputs. Furthermore, it is assumed that the valueβ/K is globally known 2 . It produces local values of σ l,j , (l, j) ∈L at each iteration, used to locally update σ as explained in Section IV. The aim of Algorithm 1 is to obtain a vector σ ∈ B |L| , such that (6) is satisfied with θ =β when σ * = σ. To achieve this, it introduces global auxiliary variables p c min , p c max ,p c , p c set ∈ R and φ ∈ B and local auxiliary variablesσ l,j ∈ B, (l, j) ∈L. Note that, as also follows from Fig. 1, (A.2), (A.7) are implemented at each load while the rest at the coordinating controller.
Intuition: Algorithm 1 aims to obtain values of p c set andp c such that the termination conditionp c (k) −β/K ≤ p c set (k) ≤ p c (k), as follows by (A.4)-(A.6), is satisfied at some finite iteration k. We show in Theorem 1 below that the termination condition (i.e. φ = 1) suffices for (6) to hold at steady state when σ * = σ. Variable p c set (k) corresponds to some price signal used to updateσ l,j (k). Its terminal value enables to obtain some ζ that satisfies (6), as explained in the proof of Theorem 1 below. The update rule forσ is intuitive, allowinĝ σ to be different than ρ only if the magnitude of its local perturbed cost per unit demand, given byc l,j /d l,j , is less than p c set . Moreover,p c (k) provides an estimate of the equilibrium value for p c, * when σ * =σ(k). The variables p c min (k) and . The latter, as already explained, enables (6) to be satisfied.
Remark 2: Note that we have opted for a distributed power command scheme, described by (5), along Algorithm 1 which includes a central coordinating controller. This choice was made for two reasons. First, the discrete nature of Algorithm 1 enables reduced bandwidth requirements in the hierarchical communication network compared to those imposed by a continuous hierarchical scheme that could potentially replace (5). Second, an interesting extension of the presented scheme would be a distributed implementation of Algorithm 1. Demonstrating that Algorithm 1 can be implemented along a distributed continuous controller for power command is a step in this direction. Extending Algorithm 1 to a fully distributed scheme introduces several challenges and is left for future work.
Remark 3: As shown in Theorem 1 below, Algorithm 1 ensures that (6) is satisfied at steady state with θ =β, which enables an optimality interpretation to be obtained. Note that (6) is based on a continuous relaxation of (4) and the corresponding KKT conditions. The design of Algorithm 1 allows to overcome the challenges associated with Problem 1 (see Section II-D).
Below, we demonstrate that Algorithm 1 terminates (i.e. φ = 1) after a finite number of iterations. Furthermore, we show that if σ * = σ * , where σ * is the value obtained for σ when Algorithm 1 terminates, then (6) is satisfied at steady state with θ =β. The latter allows an optimality interpretation of the steady state power allocation, in accordance with Proposition 1. The above are demonstrated in the following theorem, proven in the Appendix.
Remark 4: The properties of Theorem 1 could also be obtained by a centralized scheme. Such a scheme could provide demand values for on-off loads such that an ǫoptimal power allocation is enabled at steady state, using the same inputs as Algorithm 1. The advantages of the proposed hierarchic control scheme in comparison with a centralized one are multiple. Firstly, the proposed scheme does not require central knowledge of load cost coefficients and magnitudes c and d, which are assumed to be known by the loads only. Moreover, the scheme is able to accommodate changes in the load costs (i.e. different costs may be incurred at different times of the day) and the addition or removal of on-off loads without requiring to notify the system operator. By contrast, a centralized scheme would require knowledge of such changes. In addition, note that the computational effort required for Algorithm 1 to provide a close to optimal response is low since the algorithm relies on simple operations. On the other hand, a fully centralized scheme could potentially lead to a less conservative value for ǫ.

IV. CONVERGENCE ANALYSIS
In this section we provide a detailed description of the overall dynamical system as a switching system (see e.g. [35]) and use corresponding tools for its analysis. In particular, system (1), (2), (3), (5), with on-off controllable loads that switch according to Algorithm 1, can be described by the states z = (x, σ), where x = (η, ω, p M , p c , ψ) ∈ R n , n = 3|N | + |E| + |Ẽ| is the continuous state, and σ ∈ B |L| the discrete state. We also denote by Λ = R n ×B |L| the domain where the state z takes values. For convenience, we use the following compact representation to describe (1), (2), (3), (5), where f σ : R n → R n is given by (1), (2), (3), (5). Furthermore, we let t k be the time instant at which the kth iteration of Algorithm 1 is completed, satisfying t k+1 ≥ t k , k ≥ 1, assuming that the time required for each iteration of Algorithm 1 is finite. The switching state σ is given by where σ(k) is the output of the kth iteration of Algorithm 1 and σ(0) the initial value of σ. Note that, as shown in Theorem 1, Algorithm 1 terminates after a finite number of iterations, which we denote byk, and hence tk +1 is not well-defined. To resolve this, we let tk +1 = ∞.

A. Analysis of equilibria and solutions
Before investigating the stability properties of the switching system (7)-(8), we characterize its equilibria, and establish existence and completeness of solutions.
For the analysis of system (7), we will be considering its Caratheodory solutions (e.g. [35, Ch. 1]). A Caratheodory solution to (7) on an interval [0, t 1 ] is an absolutely continuous map x(t), x : [0, t 1 ] → R n that satisfies (7) for almost all t ∈ [0, t 1 ]. Caratheodory solutions are frequently used for the analysis of switching systems [35]. The use of Caratheodory solutions enables to resolve complications related to the discontinuity in the vector field as a result of the switching. For convenience, we will refer to Caratheodory solutions of (7) by just solutions. Moreover, note that a complete solution to (7) is a solution where the time domain is unbounded.
The following lemma, proven in the Appendix, establishes existence and uniqueness of complete solutions to (7)- (8).
Remark 5: Note that although we have not explicitly set a lower bound in the time between switches, chattering behavior can be excluded for (7)- (8). The latter follows from Theorem 1 which guarantees a finite number of iterations of Algorithm 1 and hence a finite number of switches for the state σ.

B. Stability analysis
In this section we provide our main convergence result for system (7)- (8). In particular, we demonstrate that solutions to (7)-(8) globally converge to the set of its equilibria.
Remark 6: Theorem 2, proven in the Appendix, provides convergence guarantees for (7)- (8), demonstrating that the stability of the power system is not compromised from the addition of the switching dynamics resulting from Algorithm 1. This, together with the absence of chattering, explained in Remark 5, demonstrates a smooth performance when controllable on-off loads are incorporated in the power network.
The following theorem, proven in the Appendix, shows that when the condition on the local droop gains from Proposition 1 holds, then solutions to (7)-(8) globally converge to an ǫoptimal point of the H-OSC problem (4).
Theorem 3: Let k j = q −1 j , j ∈ N . Then solutions to (7)-(8) globally converge to a subset of its equilibria, satisfying ω * = 0 |N | and p c, * ∈ Im(1 |N | ), that are ǫ-optimal to the H-OSC problem (4) with ǫ = 3β 2 /2K. Theorem 3 synopsizes the results of the paper, demonstrating that the implementation of Algorithm 1 does not compromise the stability of the power network and allows an optimality interpretation of the steady state power allocation when on-off loads are present. In addition, Theorem 3 ensures that the frequency converges to its nominal value at steady state, and is locally verifiable and independent of the topology of the power and communication networks. Finally, Algorithm 1 allows the addition and removal of on-off loads without requiring recalibration. Hence, all conditions of Problem 1 are satisfied.
Remark 7: The presented scheme enables decentralized stability guarantees to be provided, as demonstrated in the proof of Theorem 2. The latter enables to incorporate additional functionalities on on-off loads without compromising the stability of the power network. A notable such example, considered in [25], is the use of on-off loads for the provision of ancillary services to the power network by switching when certain frequency thresholds are exceeded. It can be shown that the presented stability properties are retained when the on-off load scheme in [25] is implemented in conjunction with Algorithm 1. This extension is omitted for brevity in presentation and to keep the focus of the paper on the optimality aspect of on-off loads.

V. SIMULATION ON THE NPCC 140-BUS SYSTEM
In this section we use the Power System Toolbox [36] to perform numerical simulations on the Northeast Power Coordinating Council (NPCC) 140-bus system, to validate our analytic results. The model used by the toolbox is more detailed than our analytic one, including voltage dynamics, line resistances, a transient reactance generator model, and high order turbine governor models 3 .
The NPCC network consists of 47 generation and 93 load buses and has a total real power of 28.55 GW. For our simulation, we considered a step increase in demand of magnitude 2 per unit (p.u. -base 100 MVA) at load buses 2 and 3 at t = 1 second. Furthermore, we considered 500 on-off loads at each of buses 1 − 20 with local load magnitudes d and local costs c selected from uniform distributions with ranges [0, 0.008]p.u. and [0, 0.1]p.u. respectively. The value of µ, used in the implementation of Algorithm 1 was also selected from a uniform distribution with range [0.005, 0.995]. In addition, we let δ = 10 −5 . The initial value of σ was selected such that (6) was satisfied for some randomly selected value for ρ ∈ B |L| . The system contains 24 buses with turbine governor generation units of which 22 were assumed to contribute in secondary frequency regulation. Their power command droop gains were selected such that the power allocation among generators and on-off controllable loads was comparable. In addition, the generator cost coefficients were selected such that q j = k −1 j at all units, in accordance with the condition in Theorem 3.
The system was tested under two different cases: (i) the system did not include any controllable loads, (ii) on-off loads implementing Algorithm 1 were included. The duration of each iteration of Algorithm 1 was assumed to be 0.3 seconds. The frequency response for these two cases at a randomly selected bus is shown on Figure 2. Figure 2 depicts a smooth frequency response for case (ii), similar to case (i), despite the additional control layer associated with Algorithm 1. Moreover, it demonstrates that the frequency converges to its nominal value, in agreement with the presented analysis.
To demonstrate the validity of the optimality analysis, the obtained results were compared with the optimal solution to (4) obtained using the Gurobi optimizer [37]. As shown in Figure 3, the proposed algorithm provides a solution that converges to the value obtained from the Gurobi optimizer, which corresponds to the global minimum of (4). In particular, it was seen that the vector σ provided from Algorithm 1 and the corresponding provided by the Gurobi optimizer were identical. The latter are in agreement with Theorem 3. By contrast, case (i) resulted in a higher cost. Note that (5) was implemented in both cases (i) and (ii). The latter allowed an optimal allocation among generating units in case (i). Implementing alternative means of generation control (e.g. integral action) in case (i) would lead to a higher cost.
To demonstrate the speed of convergence of Algorithm 1, we measured the number of iterations required for its convergence for 5000 random test cases, where at each case the values for d l,j , c l,j , (l, j) ∈L, were randomly selected from uniform distributions with ranges as described above. The results are depicted in Figure 4, which demonstrates that in the vast majority of cases (> 90%) the algorithm converges within 50 iterations. In addition, in 96% of cases the number of required iterations was less than 100, while the number of cases with more than 200 iterations amounted to less than   2%. These suggest the fast convergence of the algorithm, well within the secondary frequency control timeframe. Note that in all cases the value of µ was randomly selected from a uniform distribution with range [0.001, 0.999]. To further explore its impact on the convergence speed of the algorithm, we fixed the value of µ and compared the required number of iterations for its convergence for different values of µ on a total of 22000 random test cases. The results are demonstrated in Figure 5 which depicts the mean and 90% range (i.e. the range where 90% of the cases lied) of the number of iterations required for the convergence of Algorithm 1 for various values of µ. From Figure 5, it follows that the number of iterations decreases when µ is close to 0.5. The latter suggests that suitably choosing µ may have a significant impact to the speed of convergence of the algorithm. Note that the initialisation of the algorithm (i.e. p c min (0), p c max (0)) also affects the number of iterations which justifies why more iterations seem to be required closer to µ = 0 rather than µ = 1.
To further compare the solutions of Algorithm 1 and the Gurobi optimizer, we tested both schemes on the above described setting on 5000 cases with randomly selected parameters. The percentage cost difference in these cases is depicted on Figure 6. From Figure 6, it follows that the solutions of the two schemes coincided in more than 90% of the cases, since identical costs have been observed. Furthermore, for the remaining cases the additional cost was seen to be below the theoretical bound of ǫ, that amounted to less than 0.0001% of the total cost, confirming Theorem 3.

VI. CONCLUSION
We have considered the problem of controlling generation and on-off loads in power networks such that stability is guaranteed and a close to optimal power allocation is attained within the secondary frequency control timeframe. A mixed integer optimization problem has been formulated which ensures that the secondary frequency control objectives are satisfied at steady state. For the considered problem, analytic conditions are derived for generation and on-off loads such that an ǫoptimal allocation is enabled at equilibrium, providing a nonconservative value for ǫ. Furthermore, a hierarchical control scheme has been proposed such that the aforementioned conditions are satisfied. The combined dynamics of the physical system and the controller are jointly studied as a switching system and analytic guarantees for stability and ǫ-optimality of the equilibrium power allocation are provided. Our results are verified with numerical simulations on the NPCC 140bus network, where a large number of on-off loads has been considered and an optimality interpretation of the steady state power allocation was obtained.
Within the proof of Proposition 1, we consider a relaxed version of the H-OSC problem (4) by allowing continuous values for controllable loads. Furthermore, we relax the discrete cost functions C d l,j toĈ d l,j as follows: where (l, j) ∈L. Hence, we define the following optimization problem, called the relaxed hybrid optimal supply control problem (RH -OSC) RH -OSC: The RH-OSC problem is convex since each component of the cost function is convex. We shall make use of subgradient techniques [38,Section 23] and the KKT conditions to solve (10), as follows from Lemma 5 below, where ∂Ĉ d l,j ( σ l,j , ρ l,j ) denotes the subdifferential ofĈ d l,j at σ l,j , ρ l,j (see e.g. [38]). Lemma 5: A point (p M , σ) is a global minimum of (10) if and only if there exists λ ∈ R such that Proof of Lemma 5: The proof follows directly from applying the subgradient KKT conditions [38,Sec. 23] to (10).
Proof of Proposition 1: Let the minimum costs of the RH-OSC and H-OSC problems be denoted by C opt and C * opt respectively and note that C opt provides a lower bound to C * opt since (10) is a relaxed version of (4), allowing d c to take continuous values. Furthermore, let C * be the cost associated with (4) at some equilibrium point to (1), (2), (3), (5). It then follows that C opt ≤ C * opt ≤ C * . To prove Proposition 1, we solve the continuous optimization problem (10) using Lemma 5 and then show that when k j = q −1 j , j ∈ N , the equilibria of (1), (2), (3), (5) that satisfy (6) are ǫ-optimal to (10) which implies that they are also ǫoptimal to (4). We then show that ǫ is given by ǫ = 3θ 2 /2K.
First, note that (11c) is equivalent to where γ l,j = c l,j /d l,j . Hence, λ determines the optimum values of all on-off loads (l, j) ∈L. Furthermore, λ needs to be sufficiently large to ensure generation and demand balance, as in (11). Below, we prove Proposition 1 when λ ≥ 0, noting that the alternative case can be proven analogously.
Letting p c, * be the equilibrium power command of (7), which is equal for all buses due to (9f), it follows that when λ = p c, * , and k j = q −1 j , j ∈ N holds, then (11b) holds. Furthermore, condition (11a) follows from the summation of equilibrium equation (9g) over all j ∈ N . Hence, when λ = p c, * , if (12) is feasible, i.e. if σ l,j ∈ {0, 1}, (l, j) ∈L, then the optimal cost to (10) is equal to that of (4). Below, we quantify the additional cost incurred when (12) is not feasible or when the obtained equilibrium is different than the optimal, i.e. when C * > C * opt . Note that (11) suggests that the value of λ is uniquely determined from ℓ = 1 T |N | p L and ρ. Furthermore, note that (6) and (9) suggest that p c, * is non-increasing (respectively non-decreasing) when ζ increases (respectively decreases).
Now consider a solution (p M , σ) to (10). The optimal cost to (10), C opt , is given by where the second step follows from (11b) and the condition k j = q −1 j , j ∈ N . In analogy, the cost of (4) at an equilibrium point to (1), (2), (3), (5), C * , that satisfies (6) is given by Then, note that when (13) holds, then equilibria to (1), (2), (3), (5) that satisfy (6), also satisfy p c, * = λ +q K , (6) and (13). The above follows from the definition of q and the fact that the cost per unit demand in loads that are switched off in σ k,l and not in σ * k,l is at least ζ. Hence, the difference between the equilibrium cost and the optimal cost satisfies Simplifying (14) and using 0 ≤ λ − ζ ≤ θ/K from (13) and q ∈ [0, θ] results to C * − C opt ≤ 3θ 2 /2K and completes the proof.
Proof of Lemma 2: Summing (9g) over all j ∈ N and noting (9d) and that ω * = 0 |N | and p c, * ∈ Im(1 |N | ) from Lemma 1 suggests that First note that all references to (6) below imply its satisfaction with θ = β. A value ζ and vector σ * such that (6) and (15) are simultaneously satisfied can be obtained using the following procedure. Initially, we let ζ = 0 and obtain σ * and p c, * from (6b) and (15). If (6a) is satisfied then the proof is complete. Alternatively, if p c, * − β/K > ζ we follow the procedure below, noting that the case where ζ > p c, * can be treated analogously.
The above procedure creates a finite sequence with up to |L| iterations, where at each iteration the value of p c, * satisfying (15) decreases, noting that at each alteration the change in the magnitude of p c, * is bounded from above by β/K. The sequence is terminated since one of the following cases will be reached: (a) there exists (l, j) = arg min (k,i)∈N c k,i /d k,i , such that when σ * l,j = 0 and ζ = c l,j /d l,j then (6) is satisfied, (b) there exists (l, j) = arg min (k,i)∈N c k,i /d k,i , such that when σ * l,j = 0 then p c, * < c l,j /d l,j . In this case setting ζ = p c, * and σ * l,j = 1 and evaluating (15) enables a solution that satisfies (6), (c) σ * = 0 |L| and (6) holds with ζ = p c, * . Hence, there always exists an equilibrium to (1), (2), (3), (5) that satisfies (6).
Proof of Theorem 1: The proof is split in two parts. In part (i), we let ǫ(k) = p c max (k) − p c min (k) and show that ǫ(k + 1) ≤ νǫ(k), for some ν ∈ (0, 1). In addition, we show that there exists a finite iterationk such that the termination condition of Algorithm 1 is satisfied. In part (ii), we show that if σ * = σ * , then (6) is satisfied at steady state with θ =β.
Part (i): First note from (A.2) that allσ l,j (k), (l, j) ∈L are outputs of monotonically decreasing functions of p c set (k) for given ℓ. Thus, the term (l,j)∈L d l,jσl,j (k) is the output of a monotonically decreasing function of p c set (k). Furthermore, from (A.3),p c (k) is a monotonically increasing function of (l,j)∈L d l,jσl,j (k). Hence,p c (k) is the output of a monotonically decreasing function of p c set (k). Algorithm 1 terminates whenp c (k) −β/K ≤ p c set (k) ≤ p c (k) holds according to (A.6). When the termination condition does not hold then either: (i) p c set (k) <p c (k) −β/K and hence p c min (k) = p c set (k) from (A.4) which results to p c set (m) > p c set (k), m > k and hence top c (m) ≤p c (k), m > k, sincep c (k) is a monotonically decreasing function of p c set (k), (ii) p c set (k) >p c (k) and hence p c max (k) = p c set (k) from (A.5) which results to p c set (m) < p c set (k), m > k and hence top c (m) ≥p c (k), m > k sincep c (k) is a monotonically decreasing function of p c set (k). From the above, ǫ(k) = p c max (k)−p c min (k) trivially satisfies ǫ(k + 1) ≤ νǫ(k), where ν = max(µ, 1 − µ) < 1. The latter suggests that ǫ(k) → 0 as k → ∞ and hence the existence of some finitem where: (i) ǫ(m) < δ/2 =β − β and (ii) there exists up to one load (l, j) such thatc l,j /d l,j ∈ [p c min (m), p c max (m)]. The latter follows since the values of c l,j /d l,j are distinct. The above suggest that the difference in the values obtained for p c, * when p c set = p c min (m) and when p c set = p c max (m) is strictly less thanβ/K and hence that the termination condition, described by (A.4)-(A.6), has been satisfied at somek ≤m.
The above allow to deduce that (6) is satisfied in all cases with θ =β, which concludes the proof.
Proof of Lemma 3: The existence of an equilibrium to (7) for any σ ∈ B |L| follows directly from the equilibrium equations (9) which are linear and unbounded. The characterization of the equilibria follows in analogy to the proof of Lemma 1.
Proof of Lemma 4: First note that the trajectory of σ follows from (8) and Algorithm 1 and is independent of x. Then, consider any time interval [0, T ]. If T < t 1 then σ is uniquely given by (8a). If T ≥ t 1 there exist switching instants Q = ∪ 1≤k<m {t k } ⊂ [0, T ], m > 1 and corresponding values of σ(k), which are uniquely given from Algorithm 1. Hence σ, which relates to σ from (8), exists and is unique.
To show the existence and uniqueness of solutions for x, first note that from the global Lipschitz property of the dynamics in (1), (2), (3), (5), it follows that a unique solution exists for x for t ∈ [0, t 1 ). After the switch at t 1 the solution can be uniquely extended starting from x(t 1 ) for t ∈ [t 1 , t 2 ). The solution x(t) is unique since it results from the concatenation of unique solutions of (7). The fact that any maximal solution is complete follows from the global Lipschitz property of the vector field in (1), (2), (3), (5) and the fact that σ ∈ B |L| is bounded. Hence, a unique complete solution to (7)-(8) exists from any initial condition z(0) ∈ Λ.
Proof of Theorem 2: First note that from Theorem 1, it follows that Algorithm 1 converges after a finite amount of iterations,k, and hence there exists some finite time T = tk such that σ(t) is constant for t ≥ T . Below, we let σ * = σ(t), t ≥ T . We then consider a Lyapunov candidate function V which is demonstrated to be non-increasing within some compact set S for all t ≥ T . Then, the results in [35] allow to deduce global convergence to the set of equilibria within S, characterized by Lemma 3.
Since σ(t) is constant for t ≥ T , it can be shown thaṫ along any solution of (7). Note that V is a function of x only, and has a strict minimum at the equilibrium point x * = (η * , ω * , p M, * , ψ * , p c, * ). Hence there exists a compact set S = {(x, σ) : x ∈ Ξ and σ = σ * } for some connected neighborhood Ξ of x * , such that solutions that lie within S at some t ≥ T stay in S for all future times. The compact set Ξ, which includes x * , is given by Ξ = {(η, ω, p M , ψ, p c ) : V ≤ǭ},ǭ > 0. Note that Ξ is compact and hence from that and (17) it follows that all solutions of (7) that evolve within Ξ are bounded.
To prove Theorem 2 we make use of Lasalle's Theorem [35,Thm. 2.2]. From Lasalle's Theorem we deduce that all solutions within S at some t ≥ T converge to the largest weakly invariant subset of the set S ∩V 0 where V 0 = {z ∈ Λ :V = 0}. We now have that ifV = 0, then p M = p M, * and also ω = 0 |N | . Moreover, from (2a) we deduce that when p M j , ω j are constant at all times, then p c j must also be constant at all times. Using (9d), (9f) and ω * = 0 |N | and summing (9g) over all j ∈ N it follows that within the invariant set p c = p c, * . Furthermore, the fact that within the considered invariant set it holds that (ω, p c ) = (ω * , p c, * ) allows to deduce that within this set (η, ψ) are equal to some constant values. Hence, z = (x, σ) converges to the set of equilibrium points in S. Note that the set of equilibria of (7) within S is characterized by Lemma 3. Finally note thatǭ in Ξ can be chosen to be arbitrarily large, which allows to deduce global convergence.
Proof of Theorem 3: The proof follows directly from Proposition 1, Lemma 2, Theorem 1 and Theorem 2.