Distributed bilayered control for transient frequency safety and system stability in power grids

1 Distrib uted bilayered control for transient frequenc y safety and system stability in po wer grids Y ifu Zhang Jorge Cort ´ es Abstract —This paper considers power networks governed by swing nonlinear dynamics and subject to disturbances. W e develop a bilayered control strategy for a subset of buses that simultaneously guarantees transient frequency safety of each individual bus and asymptotic stability of the entire network. The bottom layer is a model predictive controller that, based on periodically sampled system information, optimizes control resour ces to hav e transient frequency e volve close to a safe desir ed interval. The top layer is a r eal-time controller assisting the bottom-layer controller to guarantee transient frequency safety is actually achieved. W e sho w that control signals at both layers are Lipschitz in the state and do not jeopardize stability of the network. Furthermore, we carefully characterize the information requir ements at each bus necessary to implement the controller and employ saddle-point dynamics to introduce a distributed implementation that only requir es information exchange with up to 2-hop neighbors in the power network. Simulations on the IEEE 39-bus power network illustrate our results. I . I N T R O D U C T I O N The electric power system is operated around a nominal frequency to maintain its stability and safety . Large frequency fluctuations can trigger generator relay-protection mechanisms and load shedding [2], [3], which may further jeopardize network inte grity , leading to cascading failures. Without ap- propriate operational architectures and control safeguards in place, the likelihood of such ev ents is not negligible, gi ven that the high penetration of non-rotational renew able resources provides less inertia, possibly inducing higher frequency ex- cursions [4]. These observ ations motiv ate us to dev elop control schemes to activ ely mitigate undesired transient frequency deviations under disturbances and contingencies. Specifically , we are interested in exploiting the potential benefits of dis- tributed controllers and architectures to enable plug-and-play capabilities, the efficient orchestration among the roles of the av ailable resources, and handling the coordination of lar ge numbers of them in an adaptiv e and scalable fashion. Literatur e r evie w: Power system stability is defined as the ability of regaining operating equilibrium conditions in the presence of disturbances while keeping deviations of system states within acceptable lev els [2]. A branch of research [5], [6], [7] focuses on characterizing equilibrium and conv ergence as a function of network topology , initial conditions, and sys- tem parameters, without explicitly accounting for the potential A preliminary version appeared as [1] at the 2019 American Control Conference. This work was supported by NSF A w ard CNS-1446891 and AFOSR A ward F A9550-15-1-0108. Y ifu Zhang is with The MathW orks, Inc., Natick, MA 01760, USA ( yifu.zhang19@gmail.com ). During the preparation of this work, Y ifu Zhang was affiliated with the Department of Mechanical and Aerospace Engineering, Univ ersity of California, San Diego. Jorge Cort ´ es is with the Department of Mechanical and Aerospace En- gineering, University of California, San Diego, La Jolla, CA 92093, USA ( cortes@ucsd.edu ). disruptions in power system stability caused by mechanisms that are activ ated by frequenc y excursions beyond safe limits. V arious control strategies have been proposed to improve tran- sient frequency behavior against disturbances, including iner - tial placement [8], droop coef ficient design [9], and demand- side frequency regulations [10]. Ho wever , these methods rely on some a-priori explicit frequency o vershoot estimation based on reduced-order models, and hence only provide approximate transient frequency safety guarantees. Combining the notion of control barrier [11] and L yapunov [12] functions, our previous work [13] proposes a feedback controller that meets both requirements of transient frequency safety and asymptotic sta- bility . This controller is distributed and requires no communi- cation, in the sense that each control signal regulated on an in- dividual bus only depends on neighboring system information that can be directly measured. Ho wev er , its non-optimization- based nature may cause bounded oscillations in the closed- loop system due to the lack of cooperation among control signals. Our work [14] employs a model predictive control (MPC)-based approach to address this issue, but the prediction horizon that can be used is limited by trade-offs between the discretization accuracy and the computational complexity , limiting its performance. In addition, the implementation of the MPC-based controller is only partially distributed: gi ven a set of regions in the network, a centralized controller aggregates information and determines the control actions within each region, independently of the others. Challenges in employing MPC techniques in the context of po wer networks [15], [16], [17] include the fact that, as the equilibrium point heavily depends on modeling and network parameters that cannot be precisely kno wn, it is analytically hard to establish robust sta- bilization given that the objectiv e function generally requires knowledge of the equilibrium point; the widespread use in practice of MPC with linearized models for prediction gi ven the nonlinear nature of the dynamics of po wer networks; and the processing power and information transmission, speed and reliability requirements associated with a single operator for measured state collection, online optimization, and decision making given the lar ge number of actors and v olume of data. Statement of contribution: This paper proposes a distributed controller framework implemented on buses available for con- trol that maintains network asymptotic stability and enforces transient frequency safety under disturbances. If a b us fre- quency is initially in a prescribed safe frequency interval, then it can only ev olve within the interval afterwards; otherwise, the controller leads frequency to enter the safe interval within a finite time. The proposed controller possesses a bilayer structure. The bottom layer solves periodically a finite-horizon con vex optimization problem and globally allocates control resources to minimize the ov erall control effort. The optimiza- 2 tion problem incorporates a prediction model for the system dynamics, a stability constraint, and a relaxed frequency safety constraint. The prediction model is a linearized and discretized approximation of the nonlinear continuous-time power net- work dynamics, carefully chosen to preserv e its local nature while keeping the complexity manageable. As a consequence, in the resulting con vex optimization problem, the objectiv e function can be interpreted as the sum of local control costs, and each constraint only in volv es local decision variables. This enables us to apply saddle-point dynamics to recov er its solution in a distributed fashion by allowing each bus (resp. line) to exchange system information within its neighboring buses (resp. lines). On the other hand, the top layer , as a real- time feedback controller , acts as a compensator , bridging the mismatch between the actual continuous-time po wer network dynamics and the sampled-based information employed in the bottom layer to rigorously guarantee frequency safety . The top layer control signal regulating on a generic bus only depends on physical measurements of system information within the range of its neighboring transmission lines. W e illustrate the performance of the proposed bilayered controller architecture in the IEEE 39-bus po wer network. I I . P R E L I M I N A R I E S Here we gather notation and concepts used in the paper . Notation: Let N , R , and R > , R > denote the set of nat- ural, real, nonnegati ve real, and strictly positive numbers, respectiv ely . V ariables belong to the Euclidean space unless specified otherwise. Denote by d a e as the ceiling of a ∈ R . For A ∈ R m × n , let [ A ] i and [ A ] i , j be its i th row and ( i , j ) th element, respecti vely . W e denote by A † its unique Moore- Penrose pseudoin verse and by range ( A ) its column space. For b ∈ R n , b i denotes its i th entry . Let 1 n and 0 n in R n denote the vector of all ones and zeros, respectiv ely . k · k denotes the 2- norm on R n . For any c , d ∈ N , let [ c , d ] N =  x ∈ N   c 6 x 6 d  . Denote the sign function sgn : R → {− 1 , 1 } as sgn ( a ) = 1 if a > 0, and as sgn ( a ) = − 1 if a < 0. Define the saturation function sat : R → R with limits a min < a max as sat ( a ; a max , a min ) =      a max a > a max , a min a 6 a min , a otherwise . Giv en C ⊂ R n , ∂ C denotes its boundary and C cl denotes its closure. For a point x ∈ R n and r ∈ R > , denote B r ( x ) ,  x 0 ∈ R n   k x 0 − x k 2 6 r  . Given a differentiable function l : R n → R , we let ∇ l denote its gradient. A function f : R > × R n → R n , ( t , x ) → f ( t , x ) is Lipschitz in x (uniformly in t ) if for every x 0 ∈ R n , there exist L , r > 0 such that k f ( t , x ) − f ( t , y ) k 2 6 L k x − y k 2 for any x , y ∈ B r ( x 0 ) and any t > 0. Giv en a function L : Y × Z → R , a point ( Y ∗ , Z ∗ ) ∈ Y × Z is a saddle point of L on the set Y × Z if L ( Y ∗ , Z ) 6 L ( Y ∗ , Z ∗ ) 6 L ( Y , Z ∗ ) holds for ev ery ( Y , Z ) ∈ Y × Z . For scalars a , b ∈ R , let [ a ] + b = a if b > 0, and [ a ] + b = max { a , 0 } if b 6 0. For vectors a , b ∈ R n , [ a ] + b ∈ R n is the vector whose i th component is [ a i ] + b i for e very i ∈ [ 1 , n ] N . Graph theory: W e introduce algebraic graph theory basics from [18]. An undirected graph is a pair G = ( I , E ) , where I is the vertex set and E ⊆ I × I is the edge set. A graph is connected if there exists a path between an y two vertices. W e denote by N ( i ) the set of neighbors of node i . An orientation procedure is to, for each generic edge e k ∈ E with vertices i , j , choose either i or j as the positiv e end and the other as the negativ e end. For a gi ven orientation, the incidence matrix D = ( d ki ) ∈ R m × n associated with G is defined as d ki = 1 if i is the positiv e end of e k , d ki = − 1 if i is the negati ve end of e k , and d ki = 0 otherwise. I I I . P R O B L E M S TA T E M E N T W e introduce here model for the power network and state the desired performance goals on the controller design. W e use a connected undirected graph G = ( I , E ) to rep- resent the po wer network, where I = { 1 , 2 , · · · , n } stands for the set of buses (nodes) and E = { e 1 , e 2 , · · · , e m } ⊆ I × I represents the collection of transmission lines (edges). At each bus i ∈ I , denote by ω i ∈ R , θ i ∈ R p i ∈ R , M i ∈ R > , and E i ∈ R > the shifted frequency with respect to the nominal frequency , v oltage angle, activ e power injection, inertial, and damping (droop) coef ficient, respectiv ely . Notice that we ex- plicitly allow buses to have zero inertia. W e assume at least one bus possesses strictly positi ve inertia. Giv en an arbitrary orientation procedure of G , let D ∈ R m × n be the corresponding incidence matrix. In addition, for each generic transmission line with positiv e end i and negati ve end j , denote λ i j , θ i − θ j as the voltage angle difference between node i and j ; denote b i j ∈ R > as the line susceptance. Let I u ⊂ I be the collection of bus indexes with additional control inputs. Let θ ∈ R n , ω ∈ R n , λ ∈ R m , denote the collection of θ i ’ s, ω i ’ s, and λ i j ’ s, respectiv ely . Let Y b ∈ R m × m be the diagonal matrix whose k th diagonal entry is the susceptance of the transmission line e k connecting i and j , i.e., [ Y b ] k , k = b i j . By definition, λ = D θ . (1) Let M , diag ( M 1 , M 2 , · · · , M n ) ∈ R n × n , and E , diag ( E 1 , E 2 , · · · , E n ) ∈ R n × n . The nonlinear swing dynamics of power network can be equiv alently formulated by choosing either ( θ , ω ) or ( λ , ω ) to describe the system state. Here, we use the latter one. In this case, the dynamics can be described by the follo wing differential algebraic equation [19], [20], ˙ λ ( t ) = D ω ( t ) , (2a) M ˙ ω ( t ) = − E ω ( t ) − D T Y b sin ( λ ( t )) + p + α ( t ) , (2b) where α ( t ) ∈ A ,  z ∈ R n   z w = 0 for w ∈ I \ I u  is the control signal to be designed. Furthermore, due to the trans- formation (1), one has λ ( 0 ) ∈ range ( D ) . (3) Throughout the rest of the paper, if not specified, we assume that the initial condition of system (2) satisfies (3). W e assume that the power injection p designed by the tertiary layer is balanced, i.e., 1 T n p = 0. This assumption is reasonable, given that our focus here is on the system transient frequency behavior , which instead lies within the scope of primary and secondary control. According to [5, Lemma 2], the system (2) with α ≡ 0 n has an equilibrium ( λ ∞ , 0 n ) ∈ R m + n that is locally asymptotically stable if k L † p k E , ∞ < 1 , (4) 3 where L , D T Y b D and k z k E , ∞ , max ( i , j ) ∈ E | z i − z j | for z ∈ R n . In addition, λ ∞ lies in ϒ and is unique in the closure of ϒ , where ϒ ,  λ   | λ i | < π / 2 , ∀ i ∈ [ 1 , m ] N  . Remark III.1. (Distributed dynamics). W e emphasize that the dynamics (2) is naturally distributed, i.e., the ev olution of any giv en state is fully determined by the state information from its neighbors. Specifically , for each ( i , j ) ∈ E , ˙ λ i j is determined by ω i and ω j , i.e., the states of neighbors of edge ( i , j ) ; for each i ∈ F ; ˙ ω i is determined by M i , ω i , E i , p i , α i and λ i j , b i j with ( i , j ) ∈ E that are either state, parameter , and power injections belonging to node i , or states and parameters of its neighboring edges. • Giv en a tar get subset I ω of I u , our goal is to design a distributed state-feedback controller , one per bus in I u , that maintains stability of the whole power network while cooper- ativ ely guaranteeing frequency in variance and attracti vity of nodes in I ω . Formally , the controller α should mak e the closed-loop system satisfy the following requirements: (i) F requency safety: For each i ∈ I ω , let ω i ∈ R and ¯ ω i ∈ R be lower and upper safe frequency bounds, with ω i < ¯ ω i . If ω i is initially safe, i.e., ω i ( 0 ) ∈ [ ω i , ¯ ω i ] , then we require that the entire trajectory stay within [ ω i , ¯ ω i ] . On the other hand, if ω i is initially unsafe, then we require that there exists a finite time t 0 such that ω i ( t ) ∈ [ ω i , ¯ ω i ] for e very t > t 0 . This requirement is equiv alent to asking the set [ ω i , ¯ ω i ] to be both in variant and attracti ve for each i ∈ I ω . (ii) Local asymptotic stability: The closed-loop system should preserve the asymptotic stability properties of the open- loop system (2) where α ≡ 0 n . (iii) Lipschitz continuity: The controller should be a Lips- chitz function in the state argument. This ensures the existence and uniqueness of solution for the closed-loop system and rules out discontinuities in the control signal. (iv) Economic cooperation: Each bus in I u should coop- erate with the others to reduce the ov erall control cost. (v) Distributed natur e: The controller α should be imple- mentable in distributed way , i.e., node i should be able to com- pute α i by only exchanging information with its neighboring nodes and edges. In Section IV, we introduce a centralized controller archi- tecture that meets the requirements (i)-(iv). W e later build on this architecture in Section V to provide a distributed controller that satisfies all requirements (i)-(v). I V . C E N T R A L I Z E D B I L A Y E R E D C O N T R O L L E R Here, we propose a centralized controller to address the requirements posed in Section III. Our idea for design starts from considering MPC to account for the economic coopera- tion requirement; howe ver , MPC cannot be run continuously due to the computational burden of its online optimization. W e therefore compute MPC solutions periodically . Giv en the reliance of the MPC implementation on sampled system states that are potentially outdated, we include additional compo- nents that emplo y real-time state information to tune the output of the MPC implementation and ensure stability and frequency safety . The control signal α is defined by α = α T L + α BL . (5) Roughly speaking, the bottom-layer controller α BL periodi- cally and optimally allocates control effort, while respecting a stability constraint and steering the frequency trajectories as a first step to achieve frequency in variance and attracti vity . The top-layer controller α T L , implemented in real time, slightly tunes the control trajectory generated by the bottom layer , ensuring frequency in v ariance and attracti vity . Figure 1 shows the overall structure of the closed-loop system. Interestingly , as we show later , the combination of the stability filter, low pass filter , and direct feedback control stabilizes the system regardless of what is in the MPC block. In the following, we provide detailed definitions of each of the design elements. A. Bottom-layer contr oller design W e introduce here the bottom-layer control signal α BL , which results from the combination of three components, cf. Figure 1: a MPC component, a stability filter , and a low-pass filter . The MPC component periodically samples the system state, solves an optimization problem online, and updates its output signal u M P C . The purpose of having this MPC compo- nent is to efficiently allocate control resources to achiev e the frequency safety requirement. The stability filter is designed to guarantee closed-loop asymptotic stability by enforcing monotonic decrease of an appropriate energy function (which we define later). Since ˆ u M P C is merely a piece-wise continuous signal, to av oid discontinuity in the control signal, the low- pass filter further smooths it to generate an input α BL that is continuous in time. The bottom-layer controller by itself stabilizes the system (without the need of the top layer) but does not guarantee frequency safety . This is precisely the role of the top-layer design, which based on real-time system state information, slightly tunes the control signal generated by the bottom layer to achiev e frequency safety while maintaining system stability . Note that, except for the MPC component, all other components can access real-time information. Next, we introduce each component in the bottom layer and characterize their properties. 1) MPC component: Based on the most recent sampled system information, the MPC component updates its output after solving an optimization problem online. Formally , denote { t w } w ∈ N as the collection of sampling time instants, where t w + 1 > t w > 0 holds for ev ery w ∈ N . At each sampling time t = t w , define a piece-wise continuous signal p f cst t : [ t , t + ˜ t ] → R n as the predicted value of the true power injection p for the ˜ t seconds immediately follo wing t . Note that here we particularly allow the predicted power injection to be time-varying, although its true v alue is time-in variant. For con venience of exposition, we define x , ( λ , ω , α BL ) as the augmented collection of system states (the last state comes from the lo w-pass filter component). Let x ( t w ) = ( λ ( t w ) , ω ( t w ) , α BL ( t w )) be the augmented system state value at the sampling time t w . In the predicted model, we discretize the system dynamics with time step T > 0, and denote N , d ˜ t / T e as the predicted 4 T op l a ye r c ontr ol N e two r k dy na mi c s ( 2) M P C S ta bil it y f il te r L ow - pa s s f il t e r D ir e c t f e e dba c k c ontr ol ( ( ) , ( ) , ) t t p   BL  TL  ˆ M P C u M P C u B ott om l a ye r c ontr ol Figure 1. Block diagram of the closed-loop system with the proposed controller architecture. step length. At ev ery t = t w , the MPC component solves the following optimization problem, min ˆ X , ˆ u , S g ( ˆ X , ˆ u , S ) , N ∑ k = 1  ∑ i ∈ I u c i ˆ α 2 BL , i ( k ) + ∑ i ∈ I ω d i s 2 i ( k )  s.t. F ˆ x ( k + 1 ) = A ˆ x ( k ) + B 1 ˆ p f cst ( k ) + B 2 ˆ u (6a) ˆ u ∈ A , (6b) ˆ x ( 1 ) = x ( t w ) , (6c) ω i − s i ( k ) 6 ˆ ω i ( k ) 6 ¯ ω i + s i ( k ) , ∀ i ∈ I ω , ∀ k ∈ [ 1 , N ] N , (6d) | ˆ u i | 6 ε i | α BL , i ( t w ) | , ∀ i ∈ I u . (6e) In this optimization, (6a) combines the linearized, discretized dynamics corresponding to (2) as well as the lo w-pass filter introduced later, and ˆ x , ( ˆ λ , ˆ ω , ˆ α BL ) ∈ R m + 2 n corresponds to the predicted system state. Depending on the specific dis- cretization method, one can choose different matrices F , A ∈ R ( m + 2 n ) × ( m + 2 n ) and B 1 , B 2 ∈ R ( m + 2 n ) × n (Section IV -B belo w contains a detailed discussion on discretization); ˆ p f cst t w ( k ) , p f cst t w ( t w + ( k − 1 ) T ) for ev ery k ∈ [ 1 , N ] N ; (6b) specifics the control av ailability for each bus; (6c) is the initial condi- tion; (6d) represents a soft version of the frequency safety constraint, where we penalize in the cost function the de viation of predicted frequency from its desired bounds; (6e) restricts the v alue of the control input ˆ u i ∈ R with respect to the state of the lo w-pass filter via a tunable parameter ε i > 0; finally , the objecti ve function g combines the overall cost of control effort and the penalty on the violation of the frequency safety requirement, where c i > 0 for each i ∈ I u and d i > 0 for each i ∈ I ω are design parameters. For compactness, we define ˆ X , [ ˆ x ( 1 ) , ˆ x ( 2 ) , · · · , x ( N )] , (7a) S , [ s ( 1 ) , s ( 2 ) , · · · , s ( N )] , (7b) ˆ P f cst t w , [ ˆ p f cst t w ( 1 ) , ˆ p f cst t w ( 2 ) , · · · , ˆ p f cst t w ( N )] , (7c) where for every k ∈ [ 1 , N ] N , s ( k ) is the collection of s i ( k ) ’ s ov er i ∈ I ω . W e denote by R ( G , I u , I ω , ˆ P f cst t w , x ( t w )) as the optimiza- tion problem (6) to emphasize its dependence on network topology , nodal index es with exogenous control signals, nodal indexes with transient frequency requirement, forecasted power injection, and state v alues at the sampling time. W e may simply use R if the context is clear . Also, we denote ( ˆ X ∗ , ˆ u ∗ , S ∗ ) as its optimal solution. Remark IV .1. (Selection of frequency violation penalty coef- ficient). The parameter d = { d i } i ∈ I ω in the objecti ve function plays a fundamental rule in determining how the predicted frequency can exceed the safe bounds. In the extreme case d = 0 | I ω | (i.e., no penalty for frequency violation), the MPC controller loses its functionality of adjusting frequency . As d grows, the controller ensures that the violation of the frequency safety requirement become smaller . The top-layer control introduced later adds additional input to the bottom-layer controller to ensure the frequency requirement is satisfied. • Giv en the open-loop optimization problem (6), the function u M P C corresponding to the MPC component in Figure 1 is defined as follo ws: for w ∈ N and t ∈ [ t w , t w + 1 ) , let u M P C ( t ) = ˆ u ∗ ( G , I u , I ω , ˆ P f cst t w , x ( t w )) . (8) Note the last two ar guments that ˆ u ∗ depends on: forecasted power injection v alue and state v alue of the entire network at a sampling time. T o implement (8), a straightforward idea is to ha ve one operator globally gather the abov e two values, obtain ˆ u ∗ by solving R , and finally broadcast ˆ u ∗ i to the i th node. Later in Section V, we propose an alternative distributed computation algorithm to reduce the computational burden. The next result characterizes the dependence of the controller on the sampled state values and predicted power injection. Proposition IV .2. (Piece-wise affine and continuous depen- dence of optimal solution on sampling state and pr edicted power injection). Suppose F is in vertible, then the optimiza- tion pr oblem R ( G , I u , I ω , ˆ P f cst t w , x ( t w )) in (6) has a unique optimal solution ( ˆ X ∗ , ˆ u ∗ , S ∗ ) . Furthermor e, given G , I u , and I ω , ˆ u ∗ is continuous and piece-wise af fine in ( ˆ P f cst t w , x ( t w )) , that is, there exist l ∈ N , { H ξ } l ξ = 1 , { K ξ } l ξ = 1 , { h ξ } l ξ = 1 , and { k ξ } l ξ = 1 with suitable dimensions such that ˆ u ∗ = K ξ z + k ξ , if z ∈  y | H ξ y 6 h ξ  for ξ ∈ [ 1 , l ] N (9) 5 holds for every z ∈ R ( N + 2 ) n + m , wher e z is the collection of ( ˆ P f cst t w , x ( t w )) in column-vector form. Pr oof. W e start by noting that R is feasible (hence at least one optimal solution exists) for any giv en z . This is because, gi ven a state trajectory ˆ X of (6a) with input ˆ u = 0 n and initial condi- tion (6c), choosing a sufficiently large s ( k ) for each k ∈ [ 1 , N ] N makes it satisfy constraint (6d). The uniqueness follo ws from the facts that I) g is strongly conv ex in ( ˆ u , S ) ; II) ˆ X is uniquely and linearly determined by ˆ u ; III) all constraints are linear in ( ˆ X , ˆ u , S ) . T o show continuity and piece-wise affinity , we separately consider 2 | I u | cases, depending on the sign of each { α BL , i ( t w ) } i ∈ I u . Specifically , let η , { η i } i ∈ I u ∈ { 1 , − 1 } | I u | and define B η ,  z   ( − 1 ) η i α BL , i ( t w ) > 0 , ∀ i ∈ I u  . Note that ev ery z lies in at least one of these sets and that, in any B η , the sign of each α BL , i ( t w ) with i ∈ I u is fix ed. Hence all the | I u | constraints in (6e) can be transformed into one of the following forms − ε i α BL , i ( t w ) 6 ˆ u i 6 ε i α BL , i ( t w ) if α BL , i ( t w ) > 0 , (10a) ε i α BL , i ( t w ) 6 ˆ u i 6 − ε i α BL , i ( t w ) if α BL , i ( t w ) 6 0 . (10b) Note that if α BL , i ( t w ) = 0, then ˆ u i = 0. Therefore, in ev ery B η , z appears in R in a linear fashion; hence, it is easy to re-write R into the following form: min q q T V q s.t. Gq 6 W + J η z , (11) where q is the collection of ( ˆ X , ˆ u , S ) in vector form and V  0, G , W and J η are matrices with suitable dimensions. Note that only J η depends on η . By [21, Theorem 1.12], for ev ery η ∈ {− 1 , 1 } | I u | , s ∗ is a continuous and piece- wise affine function of z whenev er z ∈ B η . Since each B η consists of only linear constraints and the union of all B η ’ s with η ∈ { 1 , − 1 } | I u | is R ( N + 2 ) n + m , one has that s ∗ is piece- wise affine in z on R ( N + 2 ) n + m . Lastly , to show the continuous dependence of s ∗ on z on R ( N + 2 ) n + m , note that since such a dependence holds on ev ery closed set B η , we only need to prove that s ∗ is unique for every z lying on the boundary shared by dif ferent B η ’ s. This holds trivially as s ∗ is unique for e very z ∈ R ( N + 2 ) n + m , which we have pro ven above. Notice that the continuity and piece-wise af finity established in Proposition IV .2 together suf fice to ensure that ˆ u ∗ is globally Lipschitz in z , and hence in the sampled system state. T o see this point, one can easily check that max ξ ∈ [ 1 , l ] N k K ξ k qualifies as a global Lipschitz constant. In addition, Proposition IV .2 also suggests an alternativ e to directly solve R without treating it as an optimization problem. Specifically , we can first compute and store { H ξ } l ξ = 1 , { K ξ } l ξ = 1 , { h ξ } l ξ = 1 , and { k ξ } l ξ = 1 , and then compute ˆ u ∗ online via (9). Howe ver , such an approach, usually called explicit MPC [22], suffers from the curse of dimensionality , in that the number of regions l gro ws exponentially fast in m + n , input size | I u | , and horizon length N . 2) Stability and low-pass filters: Here we introduce the stability and low-pass filters, explain the motiv ation behind their definitions and characterize their properties. Note that the sampling mechanism used for the MPC component ine vitably introduces delays in the bottom layer . Specifically , for any time t ∈ ( t w , t w + 1 ) , i.e., between two adjacent sampling times, u M P C ( t ) is fully determined by the old sampled system in- formation at time t w , as opposed to the current information. T o eliminate the potential negati ve effect of delay on system stability , we introduce a stability filter to enforce closed-loop stability . The low-pass filter after the stability filter simply smooths the output of the stability filter to ensure that the output of the bottom layer is continuous in time. Formally , for e very i ∈ I at any t > 0, define the stability filter as ˆ u M P C , i ( α BL ( t ) , u M P C ( t )) = sat ( u M P C , i ( t ) ; ε i | α BL , i ( t ) | , − ε i | α BL , i ( t ) | ) , (12) and define the low-pass filter as ˙ α BL , i ( t ) = − 1 τ i α BL , i ( t ) − ω i ( t ) + ˆ u M P C , i ( t ) , ∀ i ∈ I u , α BL , i ≡ 0 , ∀ i ∈ I \ I u , (13) where the tunable parameter τ i ∈ R > determines the bandwidth of the low-pass filter . In addition, although for compactness we define a stability filter for ev ery i ∈ I , one can easily see that ˆ u M P C , i ≡ 0 for ev ery i ∈ I \ I u . Both the stability and the lo w-pass filters possess a natural distributed structure: for each i ∈ I u , α BL , i only depends ω i and ˆ u M P C , i , where the latter one only depends on u M P C , i and α BL , i . This implies that to implement ˆ u M P C , i and α BL , i , it only requires local information at node i . Throughout the rest of the paper, we interchangeably use ˆ u M P C , i ( α BL ( t ) , u M P C ( t )) and ˆ u M P C , i ( t ) for simplicity . The next result establishes that ˆ u M P C is Lipschitz continuous in the system state and an important property of the bottom- layer controller α BL that we use later to establish stability . Lemma IV .3. (Lipschitz continuity and stability condition). F or the signal ˆ u M P C defined in (12) , ˆ u M P C is Lipsc hitz in system state at every sampling time t = t w with w ∈ N . Furthermor e, if α T L is Lipschitz in system state, then both α T L and α BL ar e continuous in time . Additionally , α BL , i ( t ) ˆ u M P C , i ( t ) 6 ε i α 2 BL , i ( t ) , ∀ t > 0 , ∀ i ∈ I . (14) Pr oof. If t = t w , then since | ˆ u ∗ i | 6 ε i | α BL , i ( t w ) | by (6e) and u M P C , i ( t w ) = ˆ u ∗ i for every i ∈ I u , using (12) we deduce that ˆ u M P C , i ( α BL ( t ) , u M P C ( t )) | t = t w = ˆ u ∗ i . The Lipschitz conti- nuity follows by Proposition IV .2. T o show the time-domain continuity , since ˆ u M P C is Lipschitz at every sampling point and the top-layer controller is also Lipschitz by hypothesis (we demonstrate this point later in Section IV -C), one has that the solutions of both α T L and the closed-loop system (2) exist and are unique and continuous in time. Note that u M P C in (8) is defined to be a piece-wise constant signal. One has, by (12), that ˆ u M P C is piece-wise continuous, which further makes α BL a continuous signal in time due to the lo w-pass filter . Condition (14) simply follows from the definition of saturation function. Remark IV .4. (Link between designs of the MPC component and stability filter). Note that, regardless of the MPC com- ponent output u M P C , the output of the stability filter ˆ u M P C defined in (12) always meets condition (14). This implies 6 that an y inaccuracy in the MPC component (e.g., errors in sampled state measurement, forecasted power injection, or system parameters) cannot cause instability . Howe ver , to ensure the Lipschitz continuity in Lemma IV .3, we formulate constraint (6e) employing the same coef ficient ε i in the stabil- ity filter (12). It is in this sense that both are linked. • Remark IV .5. (Continuous versus periodic sampling in the MPC component). Note that if the MPC component were to sample the system state in a continuous fashion instead, then the constraint (6e) would ensure that the output of the MPC component already satisfies the stability condition (14), and hence there would be no need for the stability filter . In this regard, the role of the stability filter is to filter out the unstable parts in u M P C caused by non-continuous sampling. • B. Discretization with sparsity pr eservation As we have introduced the dynamics of the lo w-pass and stability filters, we are now able to explicitly explain the computation of matrices F , A , B 1 and B 2 in the prediction model (6a). W e first construct a continuous-time linear model by neglecting the top-layer controller and the stability filter ( α ≈ α BL and ˆ u M P C ≈ u M P C ), and then linearizing the nonlinear dynamics in Figure 1. Our second step consists of appropri- ately discretizing this linear model. Notice that the transformation from a nonlinear continuous- time nonlinear model to a discrete one does not affect closed- loop system stability due to the presence of the stability filter . In fact, any prediction model in the MPC component cannot jeopardize stability (cf. Remark IV .4). On the other hand, such a model simplification is reasonable since α BL is designed to only slightly tune the control signal, and we hav e described in Remark IV .5 ho w the stability filter barely changes its input. W e obtain the linear model by assuming α ≈ α BL and ˆ u M P C ≈ u M P C , and approximating the dynamics in Figure 1 by ˙ λ ( t ) = D ω ( t ) , M ˙ ω ( t ) = − E ω ( t ) − D T Y b λ ( t ) + p + α BL ( t ) , M i ˙ α BL , i ( t ) = − 1 τ i α BL , i ( t ) − ω i ( t ) + ˆ u M P C , i ( t ) , ∀ i ∈ I u , α BL , i ≡ 0 , ∀ i ∈ I \ I u , (15) where the first two equations come from (2) by linearizing the nonlinear sinusoid function via sin ( Y b λ ( t )) ≈ Y b λ ( t ) . Now we re-write the abo ve linear dynamics into the compact form, ˜ G ˙ x ( t ) = ˜ Ax ( t ) + ˜ B 1 p + ˜ B 2 u M P C ( t ) , (16) for certain matrices ˜ A , ˜ B 1 , and ˜ B 2 , with ˜ A stable [20] and with ˜ G a diagonal matrix whose diagonals are 1, M i with i ∈ [ 1 , n ] N , or 0. Additionally , one can easily check that the linearized dynamics (15) and (16) preserve the locality of (2b) and (13). W e consider the following three discretization methods with step size T > 0 to construct F , A , B 1 , and B 2 matrices in (6a) approximating the continuous dynamics (16). For explanatory simplicity , we here assume ˜ G is in vertible. a) Impulse in variant discretization: F , I m + 2 n , A , e ˜ G − 1 ˜ AT , B s , Z T 0 e ˜ G − 1 ˜ A τ d τ ˜ B s , s = 1 , 2 , (17) b) Forward Euler discretization: F , ˜ G , A , T ˜ A + ˜ G , B s , T ˜ B s , s = 1 , 2 , (18) c) Backward Euler discretization: F , ˜ G − T ˜ A , A , I m + 2 n , B s , T ˜ B s , s = 1 , 2 , (19) where F should be in vertible for uniqueness of solution of the discretized dynamics. Note that with a fixed T , the impulse in v ariant and backward Euler methods usually have better approximation accuracy than the forward Euler method. In fact, since all eigen val- ues of ˜ A hav e non-positiv e real part, a basic discretization requirement is that all eigen values of F − 1 A are in the unit circle to maintain stability . One can easily prove that the impulse in variant and backward Euler discretization always meet this requirement for any T > 0, but the forward Euler method requires a sufficiently small T to preserve stability; therefore, with a same predicted time horizon ˜ t , the forward Euler method has the largest predicted step length N and hence makes the optimization problem R harder to solve. On the other hand, the backward Euler method might require a small enough T to guarantee the in vertibility of F , but numerically we hav e found this to be easily satisfiable. Therefore, we set aside the forward Euler method from our considerations of discretization. On the other hand, the impulse in variant method fails to preserve the sparsity of ˜ A , ˜ B 1 , and ˜ B 2 , which are essential for the design of distributed solvers of R . Instead, the matrices F , A , B 1 and B 2 resulting from the backward Euler discretization are all sparse. This justifies our choice, throughout the rest of the paper , of the backward Euler method for discretization. C. T op-layer contr oller design In this section we describe the top-layer controller . By design, cf. (6), the bottom-layer controller makes a trade- off between the control cost and the violation of frequency safety , and hence does not strictly guarantee the latter . This is precisely the objective of the top-layer controller: ensuring frequency safety at all times by slightly adjusting, if necessary , the effect of the bottom-layer controller . Formally , for ev ery i ∈ I ω , let ¯ γ i , γ i > 0, and ω thr i , ¯ ω thr i ∈ R with ω i < ω thr i < 0 < ¯ ω thr i < ¯ ω i . W e use the design from [13] for the top layer . For i ∈ I ω , α T L , i ( x ( t ) , p ) takes the form        min { 0 , ¯ γ i ( ¯ ω i − ω i ( t )) ω i ( t ) − ¯ ω thr i + v i ( x ( t ) , p ) } ω i ( t ) > ¯ ω thr i , 0 ω thr i 6 ω i ( t ) 6 ¯ ω thr i , max { 0 , γ i ( ω i − ω i ( t )) ω thr i − ω i ( t ) + v i ( x ( t ) , p ) } ω i ( t ) < ω thr i , (20) where v i ( x ( t ) , p ) , E i ω i ( t ) + [ D T ] i sin ( Y b λ ( t )) − p i − α BL , i ( t ) , and for i ∈ I \ I ω , simply α T L , i ≡ 0. The top-layer controller can be implemented in a decentralized fashion: for each α T L , i with i ∈ I ω on bus i , its implementation only requires the bus frequency ω i , aggregated power flow [ D T ] i sin ( Y b λ ) , po wer injection p i , and i th component of the bottom-layer signal 7 α BL , i , all of which are local to bus i . Additionally , similarly to [13], one can sho w that α T L is locally Lipschitz in x . For brevity , we may use α T L , i ( x ( t ) , p ) (respectively , v i ( x ( t ) , p ) ) and α T L , i ( t ) (respecti vely v i ( t ) ) interchangeably . Each α T L , i , with i ∈ I ω , beha ves as a passi ve and myopic transient frequency regulator without prediction capabilities. W e offer the following observations about its definition: first, α T L , i only depends on local system information and does not incorporate any global knowledge; second, α T L , i vanishes as long as the current frequency is within [ ω thr i , ¯ ω thr i ] , a subset of the safe frequency interval, with no consideration for the possibility of future large disturbances; third, α T L , i can be non- zero when the current frequency is out of [ ω thr i , ¯ ω thr i ] and hence close to the safe frequency boundaries. Howe ver , this could also lead to over -reaction, especially when ¯ γ i and γ i are small, as the disturbance may disappear suddenly , in which case even without the top-layer controller, the frequency would remain safe afterwards. As pointed out abov e, the top-layer controller only steps in if the input from the bottom-layer controller is not suf ficient to ensure frequenc y safety . D. F requency safety and local asymptotic stability Having introduced the elements of both layers in Figure 1, we are no w ready to sho w that the proposed centralized control strategy meets requirements (i)-(iv) in Section III. W e focus on the first two requirements, since we have already established the Lipschitz continuity of each individual component, and the MPC component by design takes care of the economic cooperation among the controlled buses. For the open-loop system (2) with α ≡ 0 n , under condi- tion (4), the follo wing energy function [7] is identified to prov e local asymptotic stability and estimate the region of attraction, V ( x ) , 1 2 ¯ n ∑ i = 1 M i ω 2 i + m ∑ j = 1 [ Y b ] j , j a ( λ j ) , (21) where a ( λ j ) , cos λ ∞ j − cos λ j − λ j sin λ ∞ j + λ ∞ j sin λ ∞ j for e very j ∈ [ 1 , m ] N . For notational simplicity , here we assume that the first ¯ n nodes ha ve strictly positiv e inertia, whereas the rest n − ¯ n nodes have zero inertia. Due to the extra dynamics introduced by the low-pass filter , we here consider the following energy function for the closed-loop system, ¯ V ( x ) = V ( x ) + 1 2 ∑ i 6 ¯ n , i ∈ I u M i α 2 BL , i . (22) Furthermore, define the level set T ρ ,  x   λ ∈ ϒ cl , ¯ V ( x ) 6 ρ c  , (23) where ρ > 0 and c , min ˜ λ ∈ ∂ ϒ ¯ V ( ˜ λ , 0 n , 0 n ) . Now we are ready to prov e that system (2) with the proposed controller guaran- tees frequency safety and local asymptotic stability jointly . Theorem IV .6. (Bilayer ed contr ol with stability and frequency guarantees). Under condition (4) , assume that ε i τ i < 1 for every i ∈ I u , and x ( 0 ) ∈ T , then the system (2) with the bilayer ed controller defined by (5) , (8) , (12) , (13) , and (20) satisfies (i) for any i ∈ I ω , if ω i ( 0 ) ∈ [ ω i , ¯ ω i ] , then ω i ( t ) ∈ [ ω i , ¯ ω i ] for every t > 0 ; (ii) for any i ∈ I ω , if ω i ( 0 ) / ∈ [ ω i , ¯ ω i ] , then ther e e xists t 0 such that ω i ( t ) ∈ [ ω i , ¯ ω i ] for every t > t 0 . Furthermor e, ω i ( t ) monotonically appr oaches [ ω i , ¯ ω i ] befor e entering it; (iii) if the initial state ( λ ( 0 ) , ω ( 0 ) , α BL ( 0 )) is in T ρ for some 0 < ρ < 1 , then ( λ ( t ) , ω ( t ) , α BL ( t )) stays in T ρ for all t > 0 , and con verg es to ( λ ∞ , 0 n , 0 n ) . Furthermore , α ( t ) , α T L ( t ) , α BL ( t ) , ˆ u M P C ( t ) , and u M P C ( t ) all con verg e to 0 n as t → ∞ . Pr oof. It is easy to see that statement (i) is equi valent to asking that, for an y i ∈ I ω at any t > 0, ˙ ω i ( t ) 6 0 if ω i ( t ) = ¯ ω i , (24a) ˙ ω i ( t ) > 0 if ω i ( t ) = ω i . (24b) For simplicity , we only prov e (24a), and (24b) follo ws simi- larly . Note that by (2b), (5), and (20), one has M i ˙ ω i ( t ) = − E i ω i ( t ) − [ D T ] i sin ( Y b λ ( t )) + p i + α i ( t ) = − E i ω i ( t ) − [ D T ] i sin ( Y b λ ( t )) + p i + α BL , i ( t ) + α T L , i ( t ) = − v i ( t ) + α T L , i ( t ) . (25) ω i ( t ) = ¯ ω i , then − v i ( t ) + α T L , i ( t ) = − v i ( t ) + min { 0 , v i ( t ) } 6 0; hence condition (24a) holds for ev ery i ∈ I ω with M i > 0. T o establish the result for the case when M i = 0, we reason as follows. Starting from the last line of (25), the following holds when ω i ( t ) > ¯ ω thr i , E i ω i ( t ) = E i ω i ( t ) − v i ( t ) + α T L , i ( t ) = min { E i ω i ( t ) − v i ( t ) , E i ω i ( t ) + ¯ γ i ( ¯ ω i − ω i ( t )) ω i ( t ) − ¯ ω thr i } 6 E i ω i ( t ) + ¯ γ i ( ¯ ω i − ω i ( t )) ω i ( t ) − ¯ ω thr i , and hence ¯ γ i ( ¯ ω i − ω i ( t )) ω i ( t ) − ¯ ω thr i > 0, implying that ω i ( t ) 6 ¯ ω i . Similarly , one can pro ve that ω i ( t ) > ω i for e very t > 0. Note that (ii) follo ws from (i) and (iii). This is because, for any i ∈ I , if ω i con verges to 0 ∈ ( ω i , ¯ ω i ) , there must exist a finite time t 0 such that ω i ( t 0 ) ∈ [ ω i , ¯ ω i ] , which, by (i), implies that ω i ( t ) ∈ [ ω i , ¯ ω i ] at any t > t 0 . W e then prove statement (iii), T o show the in variance of T ρ , first, it is easy to see that c > 0 by noticing that λ ∞ 6∈ ∂ ϒ , and V ( ˜ λ , 0 n , 0 n ) is non-negati ve, equaling 0 if and only if ˜ λ = λ ∞ . Next, we show that ˙ ¯ V 6 0 for e very x ∈ T ρ . W e obtain after some computations that ˙ ¯ V = − ω T ( t ) E ω ( t ) + ∑ i 6 ¯ n , i ∈ I ω ω i ( t ) α T L , i ( t ) − ∑ , i ∈ I u  1 τ i α 2 BL , i ( t ) − α BL , i ( t ) ˆ u M P C , i ( t )  . Note that by the definition of α T L in (20), ω i ( t ) α T L , i ( t ) 6 0 holds for every i ∈ I ω at e very t > 0, in that α T L , i ( t ) = 0 whenev er ω thr i 6 ω i ( t ) 6 ¯ ω thr i , and α T L , i ( t ) > 0 (reps. 6 0) if ω i ( t ) > ¯ ω thr i > 0 (respectiv ely , ω i ( t ) 6 ω thr i < 0). Therefore, together with condition (14) in Lemma IV .3, we hav e ˙ ¯ V 6 − ω T ( t ) E ω ( t ) − ∑ i ∈ I u ( 1 τ i − ε i ) α 2 BL , i ( t ) 6 0 , (26) and hence ¯ V ( x ( t )) 6 ρ c for all t > 0. Finally , by the definition of c , one can check that λ stays in ϒ cl all the time, otherwise 8 there exists some t > 0 such that λ ( t ) ∈ ∂ ϒ , resulting in ¯ V ( x ( t )) > c > ρ c . Therefore, the set T ρ is in variant. The con ver gence of state follo ws by LaSalle In variance Principle [12, Theorem 4.4]. Specifically , ω ( t ) and α BL ( t ) con- ver ge to 0 n (notice that α BL , i ≡ 0 for each i ∈ I \ I u ). Next we show that lim t → ∞ α T L , i ( t ) = 0 for ev ery i ∈ I ω , which implies that lim t → ∞ α T L ( t ) = 0 n as α T L , i ≡ 0 for each i ∈ I \ I ω . This simply follows from (20) since α T L , i ( t ) = 0 whene ver ω thr i 6 ω i ( t ) 6 ¯ ω thr i , where 0 ∈ ( ω thr i , ¯ ω thr i ) , and we have shown that lim t → ∞ ω ( t ) = 0 n . The con vergence of α ( t ) follo ws from its definition (5). Since (12) implies that | ˆ u M P C , i ( t ) | 6 ε i | α BL , i ( t ) | for ev ery i ∈ I at e very t > 0, one has lim t → ∞ ˆ u M P C ( t ) = 0 n . Finally , since ˆ u ∗ ( G , I u , I ω , ˆ P f cst t w , x ( t w )) is the optimal so- lution of the optimization problem (6), it must satisfy con- straint (6e), and since lim w → ∞ α BL ( t w ) = lim t → ∞ α BL ( t t ) = 0 n , one has lim w → ∞ ˆ u ∗ ( G , I u , I ω , ˆ P f cst t w , x ( t w )) = 0 n . Finally , the con vergence of u M P C ( t ) follo ws from its definition (8). Since the MPC component cannot jeopardize system closed- loop asymptotic stability , cf. Remark IV .4, as one can see in the proof of Theorem IV .6(iii), the con ver gence of λ ( t ) , ω ( t ) , α BL ( t ) , α ( t ) , α T L ( t ) , α BL ( t ) , and ˆ u M P C ( t ) does not require any a priori assumption on the output u M P C ( t ) of the MPC component. In the simulations, we sho w that ev en if we perturb u M P C ( t ) by intentionally shifting its output from its true v alue by a constant, the con ver gence of the remaining signals still holds. On the other hand, the conv ergence of u M P C ( t ) depends on the con ver gence of α BL ( t ) . In addition, since both ˆ u M P C ( t ) and u M P C ( t ) con verge to 0 n (and so does their dif ference), we also conclude that the stability filter ultimately lets the MPC component output signal pass, i.e., the stability filter preserves the optimality of the MPC component in the long run. One can also verify the independence between stability and the MPC component by noting that all stability results of The- orem IV .6 do not rely on any assumption on the forecasted power injection. Although the MPC component is a full-state feedback, due to this independence, Theorem IV .6 still holds if the measured state is delayed or inaccurate. This means that one could instead employ an output feedback controller by designing a state observer and feeding the estimated state into the MPC component without endangering stability . The minimal set of measured information required to realize the controller are: ω i , α BL , i , [ D T ] i sin ( Y b λ ) , and p i for ev ery i ∈ I ω . This information is used in the stability and lo w-pass filters, and the top-layer controller . Of course, inaccurate state and forecasted po wer injection lead to non-optimal control commands in the MPC component and higher cost. Remark IV .7. (F requency safety with time-varying power injection). If the power injection is time-varying, one can see from the proof of Theorem IV .6 that (i) still holds and (ii) partially holds, in the sense that the frequency would approach the safe interv al but may not enter it within a finite time. • Remark IV .8. (Independence of contr oller on equilibrium point). It should be pointed out that in Theorem IV .6, the proposed controller is able to locally stabilize the system with- out a priori knowledge on the steady-state voltage angle λ ∞ . Specifically , both α BL and α T L are not functions of λ ∞ . • Remark IV .9. (Contr ol framework without bottom layer). In our previous work [13], we have sho wn that the top-layer controller by itself makes the closed-loop system meet all requirements except for the economic cooperation. Such a lack of cooperation can be observed in two aspects. First, since α T L is only defined for nodes in I ω , those in I u \ I ω do not get in volved in controlling frequency transients. Second, the top- layer control is a non-optimization-based state feedback, where each α T L , i with i ∈ I ω is merely in charge of controlling the transient frequency for its own node i . • V . C O N T R O L L E R D E C E N T R A L I Z A T I O N The centralized bilayered controller meets the require- ments (i)-(iv) stated in Section III. In this section, we focus on the requirement (v) on the distributed implementation of the controller . While introducing each controller component in Figure 1, our discussion has shown that only the MPC compo- nent requires access to global system information, whereas all other components can be implemented in a distributed fashion. In this section, we show that by having each node and edge communicate within its 2-hop neighbors, one can solve the op- timization problem R in (6) online and hence exactly recover the MPC component ˆ u ∗ in (8). The key idea is to properly assign the decision variables in the optimization problem to each node so that the cost function can be represented as sum of local costs and the constraints can be written locally . Once this is in place, we report to saddle-point dynamics to find the solution of R in a distributed way . A. Str ong con vexification of the objective function W e start here by transforming the optimization problem R into an equivalent form whose objectiv e function is strongly con vex in all its arguments. Such property is useful later when characterizing the con ver gence properties of distributed algorithm to the optimizer . Formally , let g aug ( ˆ X , ˆ u , S ) , N − 1 ∑ k = 1 k F ˆ x ( k + 1 ) − A ˆ x ( k ) − B 1 ˆ p f cst ( k ) − B 2 ˆ u k 2 2 + N ∑ k = 1 ∑ i ∈ I c i ˆ α 2 BL , i ( k ) + ∑ i ∈ I ω d i s 2 i ( k ) ! (27) + k ˆ x ( 0 ) − x ( t w ) k 2 2 . W e denote by R aug the optimization problem with objec- tiv e function g aug and constraints giv en by (6a)-(6e). Letting Y , ( ˆ X , ˆ u , S ) ∈ R ( m + 2 n + | I ω | ) N + n , we can re-write R aug into the following compact form min Y 1 2 Y T HY + f T Y + a s.t. R 1 Y 6 r 1 , (28a) R 2 Y = r 2 , (28b) for suitable H ∈ R (( m + 2 n + | I ω | ) N + n ) × (( m + 2 n + | I ω | ) N + n ) , f ∈ R ( m + 2 n + | I ω | ) N + n , a ∈ R , R 1 ∈ R ( 2 | I ω | N + 2 | I u | ) × (( m + 2 n + | I ω | ) N + n ) , 9 R 2 ∈ R (( m + 2 n ) N + n −| I u | ) × (( m + 2 n + | I ω | ) N + n ) , r 1 ∈ R 2 | I ω | N + 2 | I u | , r 2 ∈ R ( m + 2 n ) N + n −| I u | . The next result shows the equi valence between R and R aug . Lemma V .1. (Equivalent transformation to str ong con vexity). The optimization problem R and R aug posses exactly the same optimal solution. Furthermore, if F is in vertible, then g aug is str ongly con vex in ( ˆ X , ˆ u , S ) . Pr oof. The equi valence between R and R aug follows by noting that g aug corresponds to augmenting g with equality con- straints. For notational simplicity , we assume that c i = 1 for all i ∈ I and d i = 1 for all i ∈ I ω (the proof holds for general positiv e values with minor modifications). T o show strong con vexity , one can write H as an upper-triangular block matrix, whose diagonal matrices are F T F + J T J , F T F + A T A + J T J , A T A + J T J , B T 2 B 2 , and I | I ω | N , where J ∈ R ( m + 2 n ) × n is a matrix mapping the whole state ˆ x to the partial state ˆ α BL , i.e., ˆ α BL = J ˆ x . It is easy to see that both J and B 2 are full- column-rank matrices, which, together with the inv ertibility assumption on F , implies that all fiv e matrices are positive definite. Hence, all eigen values of H are real and strictly positiv e, leading to strong con vexity of g aug , as claimed. B. Separable objective with locally expr essible constraints Next, we explain ho w the problem data defining the op- timization R aug has a structure that makes it amenable to distributed algorithmic solutions. W e start by assigning the decision variables Y = ( ˆ X , ˆ u , S ) in R aug to the nodes and edges in the netw ork. W e partition the states into voltage angle differ - ence, frequency , and low-pass filter state, i.e., ˆ x = ( ˆ λ , ˆ ω , ˆ α BL ) . For e very k ∈ [ 0 , N ] N , i ∈ [ 1 , n ] N , and j ∈ [ 1 , m ] N , we assign ω i ( k ) , ˆ u i , and ˆ α BL , i ( k ) to the i th node, and ˆ λ j ( k ) to the j th edge. F or e very i ∈ I ω , we assign s i ( k ) to the i th node. In the subsequent discussion, we say a constraint or function is local for the po wer network G if its decision variables are all from either of the following two cases: a) a node i ∈ I and its neighboring edges ( i , j ) ∈ E , and b) an edge ( i , j ) ∈ E and its neighboring nodes i and j . W e claim that (i) if F , A , B 1 and B 2 are determined by (19), then ev ery constraint in (6) is local. (ii) the objectiv e function g aug can be written as a sum of local objecti ve functions. T o see (i), note that (6b)-(6e) are a collection of constraints, each depending only on variables owned by a single node. Constraint (6a) is also local by noticing the following two points. First, the dynamics of each state in (15) is uniquely determined by the states of its neighbors. Second, we hav e shown in Section IV -B that the backward Euler discretiza- tion (19) preserv es locality . T o see (ii), first note that the sum of α 2 BL , i ( k ) (respectiv ely , s 2 i ( k ) ) over i is naturally the sum of local v ariables. Second, the two-norm square of F ˆ x ( k + 1 ) − A ˆ x ( k ) − B 1 ˆ p f cst ( k ) − B 2 ˆ u for e very k ∈ [ 1 , N − 1 ] N is the sum of square of all its m + 2 n entries, where each entry is local due to the locality of discretized dynamics. Similarly , k ˆ x ( 0 ) − x ( t w ) k 2 2 is also the sum of local variables. C. Distributed implementation via saddle-point dynamics Here we introduce a saddle-point dynamics to recov er the unique optimal solution Y ∗ of R aug in a distributed fashion. W e start from the Lagrangian of R aug L ( Y , η , µ ) = g aug ( Y ) + η T ( R 1 Y − r 1 ) + µ T ( R 2 Y − r 2 ) , (29) where η ∈ R 2 | I ω | N + 2 | I u | > 0 and µ ∈ R ( m + 2 n ) N + n −| I u | are the Lagrangian multiplier corresponding to constraints (28a) and (28b), respectively . Note that we have sho wn that a) R is feasible (cf. Proposition IV .2), b) R and R aug are equiv alent (cf. Lemma V .1), and c) all constraints in R aug are linear . These three points together imply that the refined Slater condition and strong duality hold, [23, Section 5.2.3], which further implies that at least one primal-dual solution ( Y ∗ , η ∗ , µ ∗ ) of R aug exists, and the set of primal-dual solutions is exactly the set of saddle points of L on the set R ( m + 2 n + | I ω | ) N + n × ( R 2 | I ω | N + 2 | I u | > 0 × R ( m + 2 n ) N + n −| I u | ) [23, Section 5.4.2]. There- fore, one can apply the saddle-point dynamics [24] to recover one solution ( Y ∗ , η ∗ , µ ∗ ) , where ˆ u ∗ is the MPC output signal we need. F ormally , the saddle-point dynamics of R aug is ε Z d Z d τ = − ∇ Z L ( Z , η , µ ) = − ( H Z + f + R T 1 η + R T 2 µ ) , (30a) ε η d η d τ = [ ∇ η L ( Z , η , µ )] + η = [ R 1 Z − r 1 ] + η , (30b) ε µ d µ d τ = ∇ µ L ( Z , η , µ ) = R 2 Z − r 2 , (30c) where ε Z , ε η , and ε µ are tunable positi ve scalars. Giv en the strong con ve xity of g aug , the follo wing result states the global con vergence of the dynamics (30), and its proof directly follo ws from [24, Theorem 4.2]. Theorem V .2. (Global asymptotic conver gence of saddle- point dynamics). Starting fr om any initial condition ( Z ( 0 ) , η ( 0 ) µ ( 0 )) , it holds that Z ( τ ) globally asymptotically con verg es to the unique optimal solution Y ∗ of R aug . T o conclude, we justify how the saddle-point dynamics (30) can be implemented in a distributed fashion to recover Y ∗ . W e first assign ( Z , η , µ ) to different nodes and edges. In (30), the primal v ariable Z corresponds to Y , and its assignment is exactly the same, as discussed at the be ginning of Section V -B. Since all constraints are local with respect to a node or an edge, we assign each entry of ( η , µ ) to the corresponding node or edge. With this assignment, and due to locality , the dual variables dynamics (30b) and (30c) are distributed, i.e., for each entry of η or µ , if it belongs to a node (resp., edge), then its time deriv ativ e only depends on primal and dual variables of its own and of neighboring edges (resp., nodes). On the other hand, the primal dynamics (30a) requires 2-hop communication, i.e., for each entry of Z , if it belongs to a node (resp., edge), then its time deriv ativ e depends on primal and dual v ariables of its neighboring nodes (resp., edges). Note that here we do not distinguish between communica- tion network topology and the underlying physical network topology , in that they are identical. That is to say , each node or edge needs to communicate with its neighboring nodes and edges exactly determined by the giv en physical network. In general, any communication topology that has the physical 10 topology as a subgraph will also be valid, which is a common assumption, see e.g., [25], [26], [27], [28]. It would still be possible to use an independent communication netw ork at the cost of sacrificing performance. For instance, in [29], the trade- off is to have ( m + n ) 2 agents (as opposed to ( m + n ) agents here) in total to form the communication netw ork; in [30], each agent needs to maintain an estimation of the entire optimal solution, leading to O ( mN + nN ) number of estimations in total for each agent, where N denotes the number of prediction step. Here, instead, each agent only estimates its o wn component of the optimal solution, which is of size O ( N ) . Remark V .3. (T ime scale in saddle-point dynamics). Since the MPC component updates its output at time instants { t w } w ∈ N according to (8), a requirement on the saddle-point dynam- ics (30) solving R (or equi v alently R aug ) is that it returns the optimal solution within t w + 1 − t w seconds starting from t w for ev ery w ∈ N . T o achieve this, one may tune ε Y , ε η , and ε µ to accelerate the con ver gence of the saddle-point dynamics. In practice, this corresponds to running (30) on a faster time scale, which puts requirements on the hardware regarding communication bandwidth and computation time. • Remark V .4. (Comparison with contr oller with r e gional co- or dination based on network decomposition). The proposed distributed algorithm treats each b us and transmission line as an agent, and recov ers the optimal solution by allowing each agent to exchange information only with its neighbors. In our pre vious work [1], [14], we ha ve proposed an alternati ve algorithm that does not rely on participation of e very agent at the expense of not recov ering the global optimal solution. The basic idea of this alternativ e implementation is to consider a set of regions in the network. Each region, independently of the rest, possesses its own centralized controller in charge of gathering regional information and broadcasting control signals to controllers within the region. T o account for the couplings in the dynamics, flows that connect a region and the rest of the network are assumed constant when computing the controller in each region. Although there can be nodes and edges shared by multiple regions, the control signal regulated on a shared node belongs to only one region. This implementation does not recov er the exact optimal solution and only ensures partial cooperation among the inputs. • V I . N U M E R I C A L E X A M P L E S W e verify our results on the IEEE 39-bus power network shown in Figure 2. W e run all simulations in MA TLAB 2018b in a desktop with an i7-8700k CPU@4.77GHz and 16GB DDR4 memory@3600MHz. All parameters in the po wer network dynamics (2) come from the Power System T ool- box [31]. Let I ω = { 30 , 31 , 32 , 37 } be four generator buses with transient frequency requirements. The safe frequency region is [ ω i , ¯ ω i ] = [ − 0 . 2 H z , 0 . 2 H z ] for ev ery i ∈ I ω (as ω corresponds to the shifted frequency , the safe frequency region without shifting is thus [ 59 . 8 H z , 60 . 2 H z ] ). Let { 3 , 7 , 25 } be another three non-generator buses that can provide control signals, so that I u = { 3 , 7 , 25 , 30 , 31 , 32 , 37 } . T o set up the optimization problem (6) used in the MPC component (8), we use (19) for the discretization. The controller parameters are parameter value parameter value ˜ t 10 s p f cst t ( τ ) , ∀ τ ∈ [ t , t + ˜ t ] p ( τ ) T 0 . 2 s c i , ∀ i ∈ I ω 4 N 50 c i , ∀ i ∈ I u \ I ω 1 d 100 t w , ∀ w ∈ N w ε i , ∀ i ∈ I u 1 . 9 ¯ γ i and γ i , ∀ i ∈ I ω 1 τ i , ∀ i ∈ I u 0 . 5 s ¯ ω thr i and − ω thr i 0 . 1 H z T able I C O NT RO L L E R PA RA M E T ER S . summarized in T able I. In addition, we apply the saddle-points dynamics (30) to generate the output of the MPC component in a distrib uted fashion. G8 37 37 25 25 G8 37 25 G10 30 2 2 G10 30 2 1 1 1 G1 39 39 G1 39 9 9 8 8 3 3 4 4 5 5 7 7 18 18 17 17 26 26 27 27 28 G9 29 38 38 G9 29 38 14 14 14 15 15 15 16 16 24 24 24 21 21 21 22 22 G6 35 G6 35 G7 23 36 G7 23 36 20 20 G5 34 34 G5 34 20 G5 34 19 19 G4 33 19 G4 33 G3 32 10 10 13 13 12 12 G3 32 10 13 12 6 6 11 11 G2 31 G8 37 25 G10 30 2 1 G1 39 9 8 3 4 5 7 18 17 26 27 28 G9 29 38 14 15 16 24 21 22 G6 35 G7 23 36 20 G5 34 19 G4 33 G3 32 10 13 12 6 11 G2 31 Figure 2. IEEE 39-b us power network. W e first sho w that the bilayered controller defined by (5), (8), (12), (13), (20) is able to maintain the transient frequency of selected nodes within the safe region without changing the equilibrium point (cf. Theorem IV .6(i) and (iii)). Although in the dynamics (2) we assume that the power injec- tion is constant, in simulations we perturb all non-generator nodes by a time-varying power injection. Specifically , for ev ery i ∈ [ 1 , 29 ] N , let p i ( t ) = ( 1 + δ ( t )) p i ( 0 ) where δ ( t ) =          0 . 2 sin ( π t / 50 ) if 0 6 t 6 25 , 0 . 2 if 25 < t 6 125 , 0 . 2 sin ( π ( t − 100 ) / 50 ) if 125 < t 6 150 , 0 if 150 < t . The deviation δ ( t ) p i ( 0 ) has both fast ramp-up and ramp-down periods and a long intermediate constant period. W e hav e chosen it this w ay to test the capability of the controller ag ainst both slo w-v arying and fast-v arying disturbances. Figure 3(a) shows the open-loop frequency responses of nodes 30, 31, 32, and 37 (i.e., nodes with the frequency safety requirement). All four frequency trajectories, which almost overlap with each other, exceed the lower safe frequency bound 59 . 8 H z . Howe ver , with the controller enabled, in Figure 3(b), their frequencies all ev olve within the safe region, and they all return to 60 H z as the disturbance disappears. Figure 3(c) shows the corresponding control signals. Note that, due to our specific choice of c i ’ s, the controller tends to use more non- generator control signals (i.e., α 3 , α 7 , and α 25 ) than generator ones (i.e., α 31 , α 32 , α 33 , and α 37 ). Also, note that the y split 11 into two groups and the control signals within each group possess almost the same trajectories. Next we compare the performance of the proposed con- troller with other approaches. Figures 4(a) and (b) show the frequency trajectories and control signals using the controller with regional coordination based on network decomposition proposed in [1]. As mentioned in Remark V .4, although this controller achie ves frequenc y safety , it only allows control cooperation within a limited region, instead of the entire network. This can be seen from Figure 4(b), where, with the same control cost coef ficients (cf. T able I), the two groups of control trajectories are not as uniform as those in Figure 3(c) and have a larger magnitude. Figures 4(c) and (d) are the frequency and control trajectories with only the top-layer controller , as proposed in [13], cf. Remark IV .9. Since it is a non-optimization-based control strategy , each control signal does not cooperate with others. In this specific scenario, the top-layer controller leads to fluctuations even during the time interval [25,125] s , when the disturbance is constant. This is because the top-layer controller is myopic, without further consideration for the effects of the rest of the network. The economic advantage of the proposed bilayered control can be also seen by computing the overall control cost over [0,180] s of the proposed controller, the controller in [1], and the controller in [13], which are around 163, 231 and 656, resp. Next, we e xamine the role of the bottom and top layers in determining the v alue of the input signal of our distributed controller . For node 30, Figure 5(a) shows that α BL , 30 is responsible for the larger share in the overall control signal α 30 , whereas α T L , 30 provides a slightly tuning during most of the time. If we reduce the penalty d 30 from 100 to 10, in Figure 5(b), the dominance of α BL , 30 decreases, in accordance with our discussion in Remark IV .1. On the contrary , if we raise d 30 to 1000, the contrib ution of the top layer becomes much smaller , as shown in Figure 5(c). W e further look into the bottom-layer control signals at node 30. Using the same set-up as in Figure 5(a), we plot in Figure 6(a) the MPC component output signal u M P C , 30 and the stability filter output signal ˆ u M P C , 30 . They are almost identical except for a paltry difference around 140s. Next, in Figure 6(b), we purposefully add 0.1 to u M P C , 30 , i.e., the input of the stability filter is now re-defined as u M P C , 30 + 0 . 1. Notice that ˆ u M P C , 30 , unaffected by the input shift, still conv erges to 0, which coincides with our analysis after Theorem IV .6. Fig- ure 6(c) sho ws how the saddle-point dynamics (30) con verges to the value of u M P C , 30 ( 50 ) starting from an initial guess. Here we hav e used ε Z = 5 · 10 − 4 and ε η = ε µ = 2 . 5 · 10 − 4 to ensure con vergence is attained within 1 s , cf. T able I. T o illustrate the closed-loop system performance under uncertainty , in Figure 7 we simulate three dif ferent scenar- ios. In Figure 7(a), instead of having an accurate forecasted power injection, at ev ery t > 0, we let p f cst t ( τ ) = p ( t ) for all τ ∈ [ t , t + ˜ t ] , i.e., the forecasted power injection is simply the current power injection. Note that in this case the frequencies of all four controlled nodes stay within the safe region, cf. Remark IV .7; in Figure 7(b), for each generator node (i.e., node 30 to 39), we adopt a first-order model [32] with a time constant of 5 s as the generator dynamics, and note that the frequencies still stay within the safe region most of the time; in Figure 7(c), we consider both inaccurate forecasted power injection and the generator dynamics, and the frequencies still behav e well after a short fluctuation. Lastly , we sho w that the distrib uted controller is able to steer the frequency to the safe region from unsafe initial conditions. T o do this, we consider the set-up of Figure 3 but intentionally disable the controller for the first 30 seconds. For clarity , we only show the frequency and control trajectories at node 30 in Figure 8(a). Note that the frequency quickly moves abov e the safe lo wer bound after the controller becomes activ e at t = 30s. Figure 8(b) shows the control signal, where after some brief transient, α BL , 30 still dominates the overall control signal. V I I . C O N C L U S I O N S W e have considered power networks gov erned by swing nonlinear dynamics and introduced a bilayered control strategy to regulate transient frequency in the presence of disturbances while maintaining network stability . Adopting a receding hori- zon approach, the bottom-layer controller periodically updates its output, enabling global cooperation among buses to reduce the overall control effort while respecting stability and soft frequency constraints. The top-layer controller , as a continuous state feedback controller , tunes the output of the bottom-layer control signal as required to rigorously enforces frequency safety and attractivity . W e have sho wn that the entire control structure can be implemented in a distributed fashion, where the control signal can be computed by having nodes interact with up to 2-hop neighbors in the power network. Future work will explore the optimization of the sampling sequences employed in the bottom layer to improve performance, the quantitativ e e valuation of the contrib utions of the top- and bottom-layer control signals, and the analysis of the robustness of the proposed controller against delays and saturation. R E F E R E N C E S [1] Y . Zhang and J. Cort ´ es, “Double-layered distributed transient frequency control with regional coordination, ” in American Contr ol Conference , Philadelphia, P A, Jul. 2019, pp. 658–663. [2] P . Kundur, J. Paserba, V . Ajjarapu, G. Andersson, A. Bose, C. Canizares, N. Hatziargyriou, D. Hill, A. Stankovic, C. T aylor, T . V . Cutsem, and V . V ittal, “Definition and classification of po wer system stability , ” IEEE T ransactions on P ower Systems , vol. 19, no. 2, pp. 1387–1401, 2004. [3] NERC, “Balancing and frequency control, ” North American Electric Reliability Council, T ech. Rep., 2011. [4] F . Milano, F . D ¨ orfler , G. Hug, D. J. Hill, and G. V erbi ˇ c, “Foundations and challenges of low-inertia systems, ” in P ower Systems Computation Confer ence , Dublin, Ireland, June 2018, electronic proceedings. [5] F . D ¨ orfler , M. Chertkov , and F . Bullo, “Synchronization in complex oscillator networks and smart grids, ” Pr oceedings of the National Academy of Sciences , vol. 110, no. 6, pp. 2005–2010, 2013. [6] H. D. Chiang, Direct Methods for Stability Analysis of Electric P ower Systems: Theor etical F oundation, BCU Methodologies, and Applica- tions . John Wiley and Sons, 2011. [7] T . L. V u, H. D. Nguyen, A. Megretski, J. Slotine, and K. Turitsyn, “In verse stability problem and applications to rene wables integration, ” IEEE Contr ol Systems Letters , vol. 2, no. 1, pp. 133–138, 2018. [8] J. Fang, H. Li, Y . T ang, and F . Blaabjerg, “Distributed power system virtual inertia implemented by grid-connected power conv erters, ” IEEE T ransactions on P ower Electr onics , vol. 33, no. 10, pp. 8488–8499, 2018. [9] S. S. Guggilam, C. Zhao, E. Dall’Anese, Y . C. Chen, and S. V . Dho- ple, “Optimizing DER participation in inertial and primary-frequency response, ” IEEE T ransactions on P ower Systems , vol. 33, no. 5, pp. 5194–5205, 2018. 12 (a) (b) (c) Figure 3. Frequency and control input trajectories with and without transient frequency control. Plot (a) shows the open-loop frequency responses at nodes 30, 31, 32, and 37, all exceeding the lower safe bound. The closed-loop system with the distributed control has all responses stay inside the safe region in plot (b). Plot (c) shows the corresponding control trajectories. (a) (b) (c) (d) Figure 4. Comparison of frequency and control trajectories with other approaches. Plot (a) and (b) employ the controller with regional coordination based on network decomposition proposed in [1]. Plot (c) and (d) correspond to the top-layer controller, a non-optimization-based control strategy proposed in [13]. (a) (b) (c) Figure 5. Decomposition of the control signal at node 30. Plot (a), (b), and (c) show the signals generated by the two control layers at node 30 using d 30 = 10 2 , d 30 = 10, and d 30 = 10 3 , respecti vely , as v alues for the frequency safety violation penalty coefficient in the MPC component. W ith a lar ger penalty , the bottom layer plays a more significant role in the overall control signal. (a) (b) (c) Figure 6. Decomposition of the bottom-layer control signal at node 30. Plot (a) sho ws the MPC component output signal u MPC , 30 and the stability filter output signal ˆ u MPC , 30 . Both are almost identical except for a minor discrepanc y appearing around 140s. T o make their dif ference more prominent, in plot (b), we add a constant 0.1 shift to u MPC , 30 , which does not affect the con vergence of ˆ u MPC , 30 to 0 (highlighting again the fact that the MPC component cannot jeopardize system closed- loop asymptotic stability , cf. Remark IV .4). Plot (c) shows the conv ergence of the saddle-point dynamics (30) computing u MPC , 30 ( 50 ) in 0 . 1s. [10] F . T eng, M. Aunedi, D. Pudjianto, and G. Strbac, “Benefits of demand- side response in providing frequency response service in the future GB power system, ” F r ontiers in Energy Research , vol. 3, no. 36, 2015. [11] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P . T abuada, “Control barrier functions: theory and applications, ” in Eur opean Control Conference , Naples, Italy , Jun. 2019, pp. 3420–3431. [12] H. K. Khalil, Nonlinear Systems , 3rd ed. Prentice Hall, 2002. [13] Y . Zhang and J. Cort ´ es, “Distributed transient frequency control for 13 (a) (b) (c) Figure 7. Frequency responses under inaccurate information and unmodeled dynamics. Plot (a), (b), and (c) show the frequency responses at nodes 30, 31, 32, and 37 under inaccurate forecasted power injection, first-order generator dynamics, and both, respectively . (a) (b) Figure 8. Frequency and control trajectories at node 30 when the controller is turned on after 30s. In plot (a), the frequency gradually returns to the safe region once the controller kicks in. Plot (b) shows the control signals. power networks with stability and performance guarantees, ” Automatica , vol. 105, pp. 274–285, 2019. [14] ——, “Model predictiv e control for transient frequency regulation of power networks, ” Automatica , 2020, submitted. [15] H. Jiang, J. Lin, Y . Song, and D. J. Hill, “MPC-based frequency control with demand-side participation: A case study in an isolated wind- aluminum power system, ” IEEE T ransactions on P ower Systems , vol. 30, no. 6, pp. 3327–3337, 2015. [16] A. N. V enkat, I. A. Hiskens, J. B. Rawlings, and S. J. Wright, “Distributed MPC strategies with application to power system automatic generation control, ” IEEE T ransactions on Contr ol Systems T echnology , vol. 16, no. 6, pp. 1192–1206, 2008. [17] A. Fuchs, M. Imhof, T . Demiray , and M. Morari, “Stabilization of large power systems using VSC-HVDC and model predictive control, ” IEEE T ransactions on P ower Delivery , vol. 29, no. 1, pp. 480 – 488, 2014. [18] F . Bullo, J. Cort ´ es, and S. Martinez, Distributed Control of Robotic Networks , ser . Applied Mathematics Series. Princeton University Press, 2009. [19] A. R. Bergen and D. J. Hill, “ A structure preserving model for power system stability analysis, ” IEEE Tr ansactions on P ower Apparatus and Systems , vol. 100, no. 1, pp. 25–35, 1981. [20] A. Pai, Energy Function Analysis for P ower System Stability . New Y ork: Springer, 1989. [21] F . Borrelli, Constrained Optimal Control of Linear and Hybrid Systems . New Y ork: Springer , 2003. [22] A. Alessio and B. Alberto, “ A survey on explicit model predictive control, ” in Nonlinear Model Predictive Control . Springer, 2009, pp. 345–369. [23] S. Boyd and L. V andenberghe, Conve x Optimization . Cambridge Univ ersity Press, 2004. [24] A. Cherukuri, E. Mallada, S. H. Low , and J. Cort ´ es, “The role of con vexity in saddle-point dynamics: L yapunov function and rob ustness, ” IEEE T ransactions on Automatic Control , vol. 63, no. 8, pp. 2449–2464, 2018. [25] E. Mallada, C. Zhao, and S. H. Low , “Optimal load-side control for frequency regulation in smart grids, ” IEEE Tr ansactions on Automatic Contr ol , v ol. 62, no. 12, pp. 6294–6309, 2017. [26] M. H. Nazari, Z. Costello, M. J. Feizollahi, S. Grijalva, and M. Egerst- edt, “Distributed frequency control of prosumer-based electric energy systems, ” IEEE Tr ansactions on P ower Systems , v ol. 29, pp. 2934–2942, 2014. [27] P . T rodden and A. Richards, “Cooperativ e distributed MPC of linear systems with coupled constraints, ” Automatica , vol. 49, no. 2, pp. 479– 487, 2013. [28] P . Giselsson, M. D. Doanb, T . Keviczk y , B. D. Schutter , and A. Rantzer , “ Accelerated gradient methods and dual decomposition in distributed model predictiv e control, ” Automatica , vol. 49, no. 3, pp. 829–833, 2013. [29] X. W ang, S. Mou, and B. D. O. Anderson, “Scalable, distributed algorithms for solving linear equations via double-layered networks, ” IEEE T ransactions on Automatic Control , 2020, to appear . [30] M. Zhu and S. Mart ´ ınez, “On distributed con vex optimization under inequality and equality constraints, ” IEEE T ransactions on Automatic Contr ol , v ol. 57, no. 1, pp. 151–164, 2012. [31] K. W . Cheung, J. Chow , and G. Rogers, P ower System T oolbox, v 3.0. Rensselaer Polytechnic Institute and Cherry Tree Scientific Software, 2009. [32] Z. W ang, F . Liu, S. H. Low , C. Zhao, and S. Mei, “Distributed frequency control with operational constraints, part I: Per-node po wer balance, ” IEEE T ransactions on Smart Grid , vol. 9, no. 4, pp. 1798–1811, 2018. Y ifu Zhang receiv ed the B.S. degree in automatic control from the Harbin Institute of T echnology , Heilongjiang, China, in 2014, and the Ph.D. degree in mechanical engineering from the University of California, San Diego, CA, USA, in 2019. In winter 2019, he interned at Mitsubishi Electric Research Laboratories, MA, USA. Currently he is a senior software quality engineer at The MathW orks, Inc., MA, USA. His research interests include distributed control and computation, model predictive control, adaptiv e control, data type optimization, function approximation, and neural network compression. Jorge Cort ´ es (M’02, SM’06, F’14) receiv ed the Licenciatura degree in mathematics from Univer - sidad de Zaragoza, Zaragoza, Spain, in 1997, and the Ph.D. degree in engineering mathematics from Univ ersidad Carlos III de Madrid, Madrid, Spain, in 2001. He held postdoctoral positions with the Univ ersity of T wente, T wente, 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 now a Professor in the Department of Mechanical and Aerospace Engineering, Univ ersity of California, San Diego, CA, USA. He is the author of Geometric, Control and Numerical Aspects of Nonholonomic Systems (Springer-V erlag, 2002) and co-author (together with F . Bullo and S. Mart ´ ınez) of Distributed Control of Robotic Networks (Princeton University Press, 2009). At the IEEE Control Systems Society , he has been a Distinguished Lecturer (2010-2014) and is currently its Director of Operations and an elected member (2018-2020) of its Board of Governors. His research interests include distributed control and optimization, network science, resource-aware control, nonsmooth analysis, reasoning and decision making under uncertainty , network neuroscience, and multi-agent coordination in robotic, power , and transportation networks.