Distributed Optimal Load Frequency Control Considering Nonsmooth Cost Functions

Distrib uted Optimal Load Frequency Control Considering Nonsmooth Cost Functions Zhaojian W ang a , Feng Liu a, ∗ , Changhong Zhao b , Zhiyuan Ma a , W ei W ei a a State K ey Labor atory of P ower Systems, Department of Electrical Engineering, Tsinghua University , Beijing 100084, China b National Renewable Ener gy Laboratory , Golden, CO, 80401, US Abstract This work addresses the distrib uted frequency control problem in power systems considering controllable load with a nonsmooth cost. The nonsmoothness e xists widely in power systems, such as tiered price, greatly challenging the design of distributed optimal controllers. In this regard, we first formulate an optimization problem that minimizes the nonsmooth regulation cost, where both capacity limits of controllable load and tie-line flo w are considered. Then, a distributed controller is deri v ed using the Clark generalized gradient. W e also prov e the optimality of the equilibrium of the closed-loop system as well as its asymptotic stability . Simulations carried out on the IEEE 68-bus system v erifies the e ff ecti veness of the proposed method. K e ywor ds: Nonsmooth optimization, distributed control, load frequenc y control, Clark generalized gradient. 1. Introduction W ith the proliferation of renewable generations, frequency control in po wer systems is facing a great challenge as power mismatch can fluctuate rapidly in a large amount. In this sit- uation, the con ventional centralized hierarchical control archi- tecture may not respond fast enough due to large inertia of the traditional synchronous generators [1, 2]. On the other hand, load-side controllable resources with fast response capabilities provide a ne w opportunity to frequenc y regulation [3]. In addi- tion, as controllable loads are usually dispersed geographically vast across the power system, a distributed architecture is more desirable for load frequency control than the centralized one. Recently , the so-called re verse engineering methodology is proposed by combining frequency control with optimal opera- tion problems in po wer systems [4, 5, 6, 7, 8]. Under this frame- work, distributed load frequency control is widely in vestigated [8, 9, 10, 11, 12, 13, 14]. In [8], an optimal load frequency control problem is formulated and a distrib uted controller is de- riv ed using controllable loads to realize primary frequenc y con- trol. T o eliminated the frequency deviation, the method is fur - ther extended in [9, 13] to realize a secondary load frequency control. At the same time, the tie-line power limit is consid- ered. The design approach is generalized in [11], where the spe- cific model requirement is eliminated. It only requires that the bus dynamics satisfy a passi vity condition to guarantee asymp- totic stability . In [10, 12], the operational constraints including regulation capacity limits and tie-line po wer limits are consid- ered, which guarantee both steady-state and transient capacity limit constraints. In [14], the distributed load frequency con- trol under time-varying and unknown power injection is inv es- I This work was supported by the National Natural Science F oundation of China ( No. 51677100, U1766206, No. 51621065). ∗ Corresponding author Email addr ess: lfeng@mail.tsinghua.edu.cn (Feng Liu) tigated, which can recover the nominal frequency ev en under unknown disturbances. The distributed load frequency control is of course a paid service, i.e., the system operator needs to pay for the controllable load to regulate their po wer . In the ex- isting literature, the cost of controllable load is assumed to be di ff erentiable, or equi v alently , the price of the controllable load is continuous. This is not true for a variety of cases, e.g., the price may have step changes when controllable load values are in di ff erent intervals. In such a situation, the regulation is in- herently nonsmooth, which makes existing methods di ffi cult to apply . This work designs a distributed controller for the optimal load frequency control in po wer systems, where the regulation cost function can be nonsmooth. W e relax the assumption of the objectiv e function from being di ff erentiable to nonsmooth. This work is partly motiv ated by [15]. Howe ver , di ff erent from it, we consider the interplay between the solving algorithm and the power system dynamics and prove the stability of the closed- loop system. Another di ff erence is that the objective function in [15] is strictly con ve x with respect to all decision variables. That is not necessary in our work, where some variables may not appear in the objecti ve function. In such a situation, we prov e the asymptotic con ver gence of the closed-loop system as well as the optimality of equilibrium. The rest of this paper is organized as follows. In Section II, we introduce some preliminaries and system models. Sec- tion III formulates the optimal load frequency control problem and introduces the distributed controller . In Section IV , con- ver gence of the closed-loop system and optimality of the equi- librium point are proved. W e confirm the performance of the controller via simulations on IEEE 68-b us system in Section V . Section VI concludes the paper . Pr eprint submitted to xxx September 23, 2024 2. Problem Description 2.1. Pr eliminaries and notations 2.1.1. Notations In this paper, use R n ( R n + ) to denote the n -dimensional (non- negati v e) Euclidean space. For a column vector x ∈ R n (matrix A ∈ R m × n ), x T ( A T ) denotes its transpose. For vectors x , y ∈ R n , x T y = h x , y i denotes the inner product of x , y . k x k = √ x T x denotes the Euclidean norm of x . Use 1 to denote the vector with all 1 elements. For a matrix A = [ a i j ], a i j stands for the entry in the i -th row and j -th column of A . Use Q n i = 1 Ω i to de- note the Cartesian product of the sets Ω i , i = 1 , · · · , n . Giv en a collection of y i for i in a certain set Y , y denotes the column vector y : = ( y i , i ∈ Y ) with a proper dimension, and y i as its components. 2.1.2. Pr eliminaries Let f ( x ) : R n → R be a locally Lipschitz continuous function and denote its Clarke generalized gradient by ∂ f ( x ) [16, Page 27]. For a continuous strictly conv ex function f ( x ) : R n → R , we hav e ( g x − g y ) T ( x − y ) > 0 , ∀ x , y , where g x ∈ ∂ f ( x ) and g y ∈ ∂ f ( y ). Define the projection of x onto a closed con ve x set Ω as P Ω ( x ) = arg min y ∈ Ω k x − y k (1) Use Id to denote the identity operator , i.e., Id( x ) = x , ∀ x . Define N Ω ( x ) = { v | h v , y − x i ≤ 0 , ∀ y ∈ Ω } . W e hav e P Ω ( x ) = (Id + N Ω ) − 1 ( x ) [17, Chapter 23.1]. A basic property of a projection is ( x − P Ω ( x )) T ( y − P Ω ( x )) ≤ 0 , ∀ x ∈ R n , y ∈ Ω (2) Moreov er , we also ha ve [18, Theorem 1.5.5] ( P Ω ( x ) − P Ω ( y )) T ( x − y ) ≥ k P Ω ( x ) − P Ω ( y ) k 2 (3) Define V ( x ) : = 1 2  k x − P Ω ( y ) k 2 − k x − P Ω ( x ) k 2  , and then V ( x ) is di ff erentiable and conv ex with respect to x [19, Lemma 4]. Moreover , we ha v e V ( x ) = 1 2 k P Ω ( x ) − P Ω ( y ) k 2 − ( x − P Ω ( x )) T ( P Ω ( y ) − P Ω ( x )) (4) ≥ 1 2 k P Ω ( x ) − P Ω ( y ) k 2 ≥ 0 (5) ∇ V ( x ) = P Ω ( x ) − P Ω ( y ) (6) where the inequality is due to (2). From (4), V ( x ) = 0 holds only when P Ω ( x ) = P Ω ( y ). 2.2. Network model A power network is usually composed of multiple buses, which are connected with each other through transmission lines. It can be modeled as a graph G : = ( N , E ), where N = { 0 , 1 , 2 , ... n } is the set of buses and E ⊆ N × N is the set of edges (transmission lines). Let m = |E| denote the num- ber of lines. The buses are divided into two types: generator buses, denoted by N g and load buses, denoted by N l . A gener- ator bus contains a generator ( possibly with certain aggregate load). A load bus has only load with no generator . The graph G is treated as directed with an arbitrary orientation and use ( i , j ) ∈ E or i → j interchangeably to denote a directed edge from i to j . W ithout loss of generality , we assume the graph is connected and node 0 is a reference node. The incidence matrix of the graph is denoted by C , and we have 1 T C = 0. W e adopt a second-order linearized model to describe the fre- quency dynamics of each bus. W e assume that the lines are lossless and adopt the DC po wer flo w model [10, 8]. For each bus j ∈ N , let θ j ( t ) denote the rotor angle at node j at time t and ω j ( t ) the frequency . 1 Let P l j ( t ) denote the controllable load. Let given constant P m j denote any change in power injec- tion, that occurs on the generation side or the load side, or both. Define θ i j = θ i − θ j as the angle di ff erence between bus i and j , and its compact form is denoted by θ e = ( θ i j , ( i , j ) ∈ E ). Then for each node j ∈ N , the dynamics are ˙ θ i j = ω i − ω j , j ∈ N (7a) ˙ ω j = 1 M j  P m j − P l j − D j ω j + X i : i → j B i j θ i j − X k : j → k B jk θ jk  , j ∈ N g (7b) 0 = P m j − P l j − D j ω j + X i : i → j B i j θ i j − X k : j → k B jk θ jk , j ∈ N l (7c) where M j > 0 are inertia constants, D j > 0 are damping con- stants, and B jk > 0 are line parameters that depend on the reac- tance of the line ( j , k ). The scenario is that: the system operates in a steady state at first. A certain power imbalance occurs due to variation of power injection P m j . Then controllable load accordingly changes its output to eliminate the imbalance. 3. Problem Formulation In this section, we first formulate the optimal load frequency problem with a nonsmooth objective function. Then, we pro- pose a distributed controller based on the Clark generalized gra- dient to driv e the po wer system to the optimal solution. 3.1. Optimization pr oblem The optimization problem is min P l j , φ j f ( P l ) = X j ∈N f j ( P l j ) (8a) s.t. 0 = P l j − P m j − X i : i → j B i j ( φ i − φ j ) + X k : j → k B jk ( φ j − φ k ) , j ∈ N (8b) P l j ≤ P l j ≤ P l j , j ∈ N (8c) θ i j ≤ φ i − φ j ≤ θ i j , ( i , j ) ∈ E (8d) 1 Sometimes, we also omit t for simplicity . 2 where P l j ≤ P l j are constants, denoting the lo wer and upper bound of P l j . θ i j ≤ θ i j are also constants, denoting the lower and upper bound of angle di ff erence. The first constraint is the local po wer balance. φ j is the virtual phase angle, which equals to θ j at the optimal solution. Use φ i j = φ i − φ j to denote the virtual phase angle di ff erence. In the DC power flo w , we have P i j = B i j θ i j , where P i j is the po wer of line ( i , j ). Thus, (8d) is in fact the tie-line power limit constraint. W e have the follo wing assumptions. Assumption 1. f j ( P l j ) is strictly con v ex. Assumption 2. The Slater’ s condition [20, Chapter 5.2.3] of (8) holds, i.e., problem (8) is feasible provided that the constraints are a ffi ne. Remark 1. Assumption 1 could be further relaxed, since a non- strictly con ve x function can be strictly conv exified by using a nonlinear perturbation [21]. Remark 2. Problem (8) allows the cost function f j ( P l j ) to be nonsmooth, which is required to be di ff erentiable in the existing literature [8, 9, 10, 11, 12, 13, 14]. Thus, the problem (8) is more general and suitable for a variety of real problems whose regulation costs are not smooth. A typical example is the tiered price, where the price discontinuously increases with respect to the amount of controllable load. It also should be noted that the decision variable φ j is absent in the objectiv e function of (8). It makes the paper not a trivial application of [15], i.e., the objective function is not required to be strictly conv e x to all the decision variables. It makes the conv ergence proof more challenging. Remark 3. In the existing literature, the controller usually in- volv es the projection of a gradient onto a con v ex set. If the objectiv e function is nonsmooth, it becomes the projection of a subdi ff erential set onto a con v ex set. In this situation, the exis- tence of trajectories is not guaranteed [15, 22], which makes ex- isting load frequenc y control methods inapplicable to the nons- mooth case. 3.2. Contr oller Design T o help the controller design, we make a modification on the problem (8). min P l j , φ j f ( P l ) = X j ∈N f j ( P l j ) + 1 2 X j ∈N z 2 j (9a) s.t. (8b) , (8c) , (8d) (9b) where z j = P l j − P m j − P i : i → j B i j φ i j + P k : j → k B jk φ jk . For any feasible solution to (8), z j = 0. Thus, (8) and (9) hav e same solutions. Define the sets Ω j : =  P l j | P l j ≤ P l j ≤ P l j  , Ω = Y n j = 1 Ω j (10) Then, we giv e the controller for each controllable load, which is denoted by OLC. ˙ d j ∈ n p : p = − d j + P l j + ω j − g j ( P l j ) − z j − µ j , g j ( P l j ) ∈ ∂ f j ( P l j ) o (11a) ˙ µ j = P l j − P m j − X i : i → j B i j φ i j + X k : j → k B jk φ jk (11b) ˙ φ j = X i : i → j B i j ( µ i − µ j ) − X k : j → k B jk ( µ j − µ k ) − X ( i , j ) ∈E η − i j + X ( j , k ) ∈E η − jk + X ( i , j ) ∈E η + i j − X ( j , k ) ∈E η + jk + X i : i → j B i j ( z i − z j ) − X k : j → k B jk ( z j − z k ) (11c) ˙ ϕ + i j = − ϕ + i j + η + i j + φ i j − θ i j (11d) ˙ ϕ − i j = − ϕ − i j + η − i j + θ i j − φ i j (11e) P l j = P Ω j  d j  (11f ) η + i j = P R +  ϕ + i j  (11g) η − i j = P R +  ϕ − i j  (11h) Combining with the power system dynamics, we have the closed-loop system (7), (11). Remark 4 (Load demand estimate) . In power systems, the load demand P l j is di ffi cult to measure. Similar to [9, 13, 12], P l j − P m j in (11b) can be substituted equiv alently in following ways. For j ∈ N g , P l j − P m j = − M j ˙ ω j − D j ω j + X i : i → j P i j − X k : j → k P jk For j ∈ N l , P l j − P m j = − D j ω j + X i : i → j P i j − X k : j → k P jk In this way , the measurement of load demand P l j is a voided. W e only need to measure ω j , P i j , which are much easier to re- alize. Moreover , the po wer loss can be treated as unknown load demand, which can be also considered by this method. 4. Optimality and Con vergence In this section, we address the optimality of the equilibrium point and the con v ergence of the closed-loop system. 4.1. Optimality Denote x = ( θ, ω g , d , µ, φ, ϕ − , ϕ + ) and y = ( x , P l , η + , η − ). Let x ∗ = ( θ ∗ , ω ∗ g , d ∗ , µ ∗ , φ ∗ , ϕ −∗ , ϕ + ∗ ) be an equilibrium of the closed-loop system (7), (11). Then, there exists g ( P l ∗ ) ∈ ∂ f ( P l ∗ ) such that 0 = C T ω ∗ (12a) 0 = P m − P l ∗ − D ω ∗ − C B C T θ ∗ (12b) 0 = − d ∗ + P l ∗ + ω ∗ − g ( P l ∗ ) − µ ∗ (12c) 0 = P l ∗ − P m + C B C T φ ∗ (12d) 0 = − C B C T µ ∗ − C η −∗ + C η + ∗ (12e) 0 = − ϕ + ∗ + η + ∗ − C T φ ∗ − θ (12f ) 0 = − ϕ −∗ + η −∗ + θ + C T φ ∗ (12g) P l ∗ = P Ω ( d ∗ ) (12h) 3 η + ∗ = P R m +  ϕ + ∗  (12i) η −∗ = P R m +  ϕ −∗  (12j) Now , we introduce the properties of the equilibrium points. Theorem 1. Suppose Assumptions 1 and 2 hold. W e have 1. The nominal frequency is restored, i.e., ω ∗ j = 0 for all j ∈ N . 2. If x ∗ is an equilibrium point of (7), (11), then ( P l ∗ , φ ∗ ) is an optimal solution to (8) and ( µ ∗ , η + ∗ , η −∗ ) is an optimal solution to its dual problem. 3. φ ∗ i j = θ ∗ i j for all ( i , j ) ∈ E . Moreover , the line limits are satisfied by x ∗ , implying θ i j ≤ θ ∗ i j ≤ θ i j on ev ery tie line ( i , j ) ∈ E . 4. At the equilibrium, ( θ ∗ , φ ∗ , ω ∗ g , P l ∗ ) is unique, with ( θ ∗ , φ ∗ ) being unique up to (equilibrium) reference angles ( θ 0 , φ 0 ). Pr oof. 1) From (12b) and (12d), we hav e 1 T D ω ∗ = 0. From (12a), we have ω ∗ = ω 0 · 1 with a constant ω 0 . As D is a diagonal positiv e definite matrix, we ha ve ω 0 = 0. 2) From (12c) and (12f)-(12j), we hav e P l ∗ = P Ω  P l ∗ − g ( P l ∗ ) − µ ∗  (13a) η + ∗ = P R m +  η + ∗ − C T φ ∗ − θ  (13b) η −∗ = P R m +  η −∗ + θ + C T φ ∗  (13c) or equiv alently , − g ( P l ∗ ) − µ ∗ ∈ N Ω ( P l ∗ ) (14a) − C T φ ∗ − θ ∈ N R m + ( η + ∗ ) (14b) θ + C T φ ∗ ∈ N R m + ( η −∗ ) (14c) By the KKT condition in [23, Theorem 3.34], (12d), (12e) and (14) coincide with the KKT optimality condition of the problem (8). Then, we have this assertion. 3) From (12b) and (12d), we have C BC T ( θ ∗ − φ ∗ ) = 0, which holds for any incidence matrix C . Thus, we have θ ∗ − φ ∗ = c 0 · 1 with a constant c 0 . Then, we have θ ∗ i j − φ ∗ i j = 0. Moreover , by 2), we know θ i j ≤ φ ∗ i j ≤ θ i j , which implies that θ i j ≤ θ ∗ i j ≤ θ i j . 4) P l ∗ is unique because the objectiv e function in (8a) is strictly con ve x in P l . ω ∗ is unique due to ω ∗ = 0. By (12d), we know φ ∗ is unique modulo a rigid (uniform) rotation of all angles. Since θ ∗ − φ ∗ = c 0 · 1 , it implies that θ ∗ is also unique modulo a rigid rotation. This proves the uniqueness of ( θ ∗ , φ ∗ , ω ∗ g , P l ∗ ). 4.2. Con ver g ence Define the function V ( x ) = V 1 ( x ) + V 2 ( x ) (15) where V 1 ( x ) = 1 2    P l − P l ∗    2 + 1 2 k µ − µ ∗ k 2 + 1 2 k φ − φ ∗ k 2 + 1 2    η + − η + ∗    2 + 1 2  θ e − θ ∗ e  T B  θ e − θ ∗ e  + 1 2    η − − η −∗    2 + 1 2  ω g − ω ∗ g  T M  ω g − ω ∗ g  (16) V 2 ( x ) = − ( d − P l ) T ( P l ∗ − P l ) − ( ϕ + − η + ) T ( η + ∗ − η + ) − ( ϕ − − η − ) T ( η −∗ − η − ) (17) Then, we hav e the follo wing result about V ( x ). Lemma 2. Suppose Assumptions 1 and 2 hold. Then the func- tion V ( x ) has following properties 1. V ( x ) ≥ 0 and V ( x ) = 0 holds only at the equilibrium point. 2. The time deriv ati v e of V ( x ( t )) satisfies ˙ V ( x ( t )) ≤ 0. Pr oof. 1) By (2), we kno w that V 2 ( x ) ≥ 0. From (15), (16) and (17), we know V ( x ) ≥ 0 and V ( x ) = 0 holds only at the equilibrium point. 2) By (6), the gradient of V is ∇ V =                               ∇ d V ∇ µ V ∇ φ V ∇ η + V ∇ η − V ∇ θ V ∇ ω g V                               =                                P l − P l ∗ µ − µ ∗ φ − φ ∗ η + − η + ∗ η − − η −∗ B  θ e − θ ∗ e  M  ω g − ω ∗ g                                 (18) Then, there is g ( P l ) ∈ ∂ f ( P l ) such that the time deriv ativ e of V is ˙ V = ( P l − P l ∗ ) T ( − d + P l + ω − g ( P l ) − z − µ ) + ( µ − µ ∗ ) T ( P l − P m + C B C T φ ) +  θ e − θ ∗ e  T B C T ω + ( φ − φ ∗ ) T  − C BC T µ − C η − + C η + − C B C T z  +  η + − η + ∗  T  − ϕ + + η + − C T φ − θ  +  η − − η −∗  T  − ϕ − + η − + θ + C T φ  + ( ω − ω ∗ ) T ( P m − P l − D ω − C B θ e ) (19) where the last item is due to the fact that, for each j ∈ N l 0 =  ω j − ω ∗ j   P m j − P l j − D j ω j + X i : i → j B i j θ i j − X k : j → k B jk θ jk  (20) Combing (19) and (12), we hav e ˙ V = ( ˜ P l ) T ( − ˜ d + ˜ P l + ˜ ω − g ( P l ) + g ( P l ∗ ) − ˜ z − ˜ µ ) + ˜ φ T  − C BC T ˜ µ − C ˜ η − + C ˜ η + − C B C T ˜ z  + ˜ µ T ( ˜ P l + C B C T ˜ φ ) +  ˜ η +  T  − ˜ ϕ + + ˜ η + − C T ˜ φ  +  ˜ η −  T  − ˜ ϕ − + ˜ η − + C T ˜ φ  + ˜ θ T B C T ˜ ω + ˜ ω T ( − ˜ P l − D ˜ ω − C B ˜ θ ) = − ( P l − P l ∗ ) T ( d − d ∗ ) +    P l − P l ∗    2 (21a) − ( η + − η + ∗ ) T ( ϕ + − ϕ + ∗ ) +    η + − η + ∗    2 (21b) − ( η − − η −∗ ) T ( ϕ − − ϕ −∗ ) +    η − − η −∗    2 (21c) − ( P l − P l ∗ ) T  g ( P l ) − g ( P l ∗ )  − ˜ ω T D ˜ ω (21d) − ( ˜ P l ) T ˜ z − ˜ φ T C B C T ˜ z (21e) 4 where ˜ x = x − x ∗ . By the [18, Theorem 1.5.5], we hav e − ( P l − P l ∗ ) T ( d − d ∗ ) +    P l − P l ∗    2 ≤ 0 (22) Similarly , − ( η + − η + ∗ ) T ( ϕ + − ϕ + ∗ ) + k η + − η + ∗ k 2 ≤ 0 and − ( η − − η −∗ ) T ( ϕ − − ϕ −∗ ) + k η − − η −∗ k 2 ≤ 0 also hold. The con v exity of f implies that − ( P l − P l ∗ ) T  g ( P l ) − g ( P l ∗ )  ≤ 0 (23) W e also have − ˜ ω T D ˜ ω ≤ 0 because D is positiv e definite. In addition, − ( ˜ P l ) T ˜ z − ˜ φ T C B C T ˜ z = − ˜ z T · ˜ z ≤ 0 (24) Then, (21a)-(21e) are all nonpositiv e, i.e., ˙ V ( x ( t )) ≤ 0. The following result shows the stability of the closed-loop system (7), (11). Theorem 3. Suppose Assumptions 1 and 2 hold. Then the tra- jectory of the closed loop system (7), (11) has following prop- erties 1. ( x ( t ) , P l ( t ) , η + ( t ) , η − ( t )) is bounded. 2. ( x ( t ) , P l ( t ) , η + ( t ) , η − ( t )) con v erges to equilibrium of the closed-loop system (7), (11). 3. The con ver gence of x ( t ) is a point, i.e., x ( t ) → x ∗ as t → ∞ for some equilibrium point x ∗ . Pr oof. 1) From Lemma 2, we kno w that ( θ ( t ) , ω g ( t ) , µ ( t ) , φ ( t ) , P l ( t ) , η + ( t ) , η − ( t )) is bounded. By (7c), ω l ( t ) is also bounded. Since ∂ f ( P l ) is compact, there exists a constant a 1 such that    P l + ω − g ( P l ) − µ    < a 1 (25) Define following function ˜ V d ( d ) = 1 2 k d k 2 (26) The time deriv ati v e of ˜ V d ( d ) along the closed-loop system is ˙ ˜ V d = d T ( − d + P l + ω − g ( P l ) − µ ) = − k d k 2 + d T ( P l + ω − g ( P l ) − µ ) ≤ − k d k 2 + a 1 k d k = − 2 ˜ V d + a 1 q 2 ˜ V d (27) Thus, ˜ V d ( d ( t )) , t ≥ 0 is bounded, so is d ( t ) , t ≥ 0. Similarly , we can also hav e that ϕ + ( t ) , t ≥ 0 and ϕ − ( t ) , t ≥ 0 are bounded. 2) By the inv ariance principle in [24, Theorem 2], we know that the trajectory x ( t ) con verges to the largest weakly in variant subset W ∗ contained in W : = { x | ˙ V ( x ) = 0 } , i.e., once a trajectory enters this subset, it will nev er departure from it. From ˜ ω T D ˜ ω = 0, we know ˜ ω = 0, i.e., ω ( t ) = ω ∗ . Note that ( P l − P l ∗ ) T  g ( P l ) − g ( P l ∗ )  > 0 if x , x ∗ due to the strict con v exity of f ( P l ). Thus, we hav e P l ( t ) = P l ∗ in the set W ∗ . Moreov er , from (24), we have ˜ P l ( t ) = C BC T ˜ φ ( t ), or 0 = ˙ ˜ P l ( t ) = C B C T ˙ ˜ φ ( t ), which implies that ˙ ˜ φ ( t ) = 0. Then, we have ˙ µ j ( t ) = 0 from (11b). Similarly , ˙ d j ( t ) = 0 by (11a). Up to now , we kno w x ( t ) is constant except ϕ − ( t ) , ϕ + ( t ) for t → ∞ . Moreov er , the equality in (3) holds only when P Ω ( x ) = P Ω ( y ) or x = P Ω ( x ) and y = P Ω ( y ). Thus, for η + ( t ) , η − ( t ) , t → ∞ , there are four combinations: 1. η + ( t ) = η + ∗ and η − ( t ) = η −∗ ; 2. η + ( t ) = η + ∗ and ϕ − ( t ) = P R m + ( ϕ − ( t ) ) = η − ( t ), ϕ −∗ = P R m + ( ϕ −∗ ) = η −∗ ; 3. ϕ + ( t ) = P R m + ( ϕ + ( t ) ) = η + ( t ) , ϕ + ∗ = P R m + ( ϕ + ∗ ) = η + ∗ and η − = η −∗ ; 4. ϕ + ( t ) = P R m + ( ϕ + ( t ) ) = η + ( t ) , ϕ + ∗ = P R m + ( ϕ + ∗ ) = η + ∗ and ϕ − ( t ) = P R m + ( ϕ − ( t ) ) = η − ( t ) , ϕ −∗ = P R m + ( ϕ −∗ ) = η −∗ . Thus, ( x ( t ) , P l ( t ) , η + ( t ) , η − ( t )) con v erges to equilibrium of the closed-loop system. 3) Fix any initial state x (0) and consider the trajectory ( x ( t ) , t ≥ 0) of the closed-loop system. As x ( t ) is bounded, there exists an infinite sequence of time instants t k such that x ( t k ) → ˆ x ∗ as t k → ∞ , for some ˆ x ∗ ∈ W ∗ . Using this specific equlibrium point ˆ x ∗ in the definition of V , we hav e V ∗ = lim t →∞ V ( x ( t )) = lim t k →∞ V ( x ( t k )) = lim x ( t k ) → ˆ x ∗ V 2  x ( t k )  = V 2 ( ˆ x ∗ ) = 0 Here, the first equality uses the fact that V ( t ) is nonincreasing in t while lo wer-bounded, and therefore must render a limit v alue V ∗ ; the second equality uses the fact that t k is the infinite subse- quence of t ; the third equality uses the fact that x ( t ) is absolutely continuous in t ; the fourth equality is due to the continuity of V ( x ), and the last equality holds as ˆ x ∗ is an equilibrium point of V . The quadratic part V 1 implies that ( θ , ω g , P l , µ, φ, η − , η + ) → ( θ ∗ , ω ∗ g , P l ∗ , µ ∗ , φ ∗ , η −∗ , η + ∗ ) as t → ∞ . Moreover , from (12c), (12f) and (12g), we can get the corresponding d ∗ , ϕ −∗ , ϕ + ∗ . This completes the proof. 5. Case studies 5.1. T est system In this section, the IEEE 68-bus New England / New Y ork in- terconnection test system [8] is utilized to illustrate the perfor- mance of the proposed controller . The diagram of the 68-bus system is gi ven in Fig.1. W e run the simulation on Matlab us- ing the Po wer System T oolbox [25]. Although the linear model is used in the analysis, the simulation model is much more de- tailed and realistic. The generator includes a two-axis subtran- sient reactance model, IEEE type DC1 exciter model, and a classical power system stabilizer model. AC (nonlinear) po wer flows are utilized, including non-zero line resistances. The up- per bound of P l j is the load demand value at each bus. Detailed simulation model including parameter values can be found in the data files of the toolbox. The objectiv e function of each controllable load is f ( P l j ) =                 P l j  2 − 0 . 02 , P l j ≤ − 0 . 2 1 2  P l j  2 , − 0 . 2 < P l j ≤ 0 . 2  P l j  2 − 0 . 02 , 0 . 2 < P l j (28) It can be verified that f ( P l j ) is continuous, strictly con ve x and nonsmooth. 5 1 2 25 26 28 29 G9 61 G8 60 G1 53 47 48 40 41 G14 66 30 31 G10 62 42 G15 67 46 38 49 52 G16 68 50 51 33 32 G11 63 34 35 45 44 39 36 G12 64 43 37 G13 65 9 8 7 6 G2 54 4 5 3 18 17 27 16 15 14 12 11 G3 55 13 10 24 19 G5 57 20 G4 56 21 G6 58 22 G7 59 23 Figure 1: IEEE 68-bus system 59.7 59.75 59.8 59.85 59.9 59.95 60 60.05 60.1 Frequency (Hz) AGC OLC 0 20 40 60 80 100 120 140 160 180 Time (s) Figure 2: Frequency dynamics under A GC and OLC 5.2. Simulation results W e consider the follo wing scenario: at t = 1s, there is a step change of (3 . 5 , 3 . 5 , 3 . 5 , 3 . 5 , 3 . 5 , 7)p.u. load demand at buses 4, 8, 20, 37, 42, and 52 respectively . Neither the original load demand nor its change is known. The load estimate method in Remark 4 is utilized. At first, we do not set limits to the tie-line po wer . In this subsection, we analyze the dynamic performance of the closed- loop system under the proposed controller OLC. In addition, automatic generation control (A GC) is tested in the same sce- nario as a benchmark. The setting of A GC is the same as that in [13]. The frequency dynamics under OLC and A GC are giv en in Fig.2. It is shown that both A GC and OLC can recover the frequency to the nominal v alue. They also hav e similar fre- quency nadir . Compared with A GC, the frequency under OLC has faster con v er gence speed. The dynamics of µ and controllable load are illustrated in Fig.3. If tie-line power is not considered, then η −∗ = η + ∗ = 0. From (12e), we know that C BC T µ ∗ = 0, i.e., µ of each bus will con v erge to the same value. This is sho wn in the left of Fig.3, where µ of buses 1 ∼ 5 is given. The right part of Fig.3 illus- trates the dynamics of controllable load at these buses. They also con ver ge to the same value except bus 2 and bus 5 as the objectiv e function is identical. As there is no controllable load on buses 2 and 5, their value is always zero. This validates the correctness of the proposed method. W e further take the tie-line limits into account. In this sce- nario, the acti v e power limit of line (1 , 2) is set to be 0. Then, its tie-line power dynamics with and without limit is given in Fig.4. The activ e po wer is − 2 . 1p.u. if there is no limit. It decreases to -0.5 0 0.5 1 1.5 2 Bus1 Bus3 Bus4 Bus7 Bus8 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 Controllable load (p.u.)  0 2 0 4 0 6 0 8 0 100 120 140 160 180 Time (s) 0 2 0 4 0 6 0 8 0 100 120 140 160 180 Time (s) Bus1 Bus3 Bus4 Bus7 Bus8 Bus1 Bus2 Bus3 Bus4 Bus5 Bus1 Bus2 Bus3 Bus4 Bus5 Controllable load (p.u.) -1 -0.8 -0.6 -0.4 -0.2 0 0.2 -0.5 0 0.5 1 1.5 2 0 20 40 60 80 100 120 140 160 180 Time (s) 0 20 40 60 80 100 120 140 160 180 Time (s)  Figure 3: Dynamics of µ and controllable load 0 20 40 60 80 100 120 140 160 180 -4 -3 -2 -1 0 1 2 3 4 5 6 Time (s) Active power of line 1-30 (p.u.) Without line control With line control -5 -4 -3 -2 -1 0 1 2 3 4 5 Active power of line (1,2) (p.u.) Without line limit With line limit 0 20 40 60 80 100 120 140 160 180 Time (s) Figure 4: Active po wer dynamics of line 1 and controllable load zero if the limit is considered, which verifies the e ff ectiv eness of the line power control. If some tie-line limit is reached, the µ of buses near from this line will diver ge. This is illustrated in the left part of Fig.5, where µ 1 , µ 3 , µ 4 are all di ff erent. W e also find that the µ of buses far from the line still conv erges to the same value, which is given in the right part of Fig.5. It is sho wn that µ 19 ∼ µ 23 con v erge to 0 . 83. 6. Conclusion In this paper , we in vestig ate the distributed load frequency control in power systems when the regulation cost is nons- mooth. In our formulation, both capacity limits of controllable load and line flow are considered. A distributed controller is designed, where the Clark generalized gradient is utilized to ad- dress the nonsmoothness of the objective function. In addition, we prove the optimality of the equilibrium of the closed-loop system as well as its asymptotic stability . Moreover , it is also prov ed that the conv er gence is to a specific point. Finally , nu- merical experiments on the IEEE 68-bus system show that the frequency is recovered to the nominal value. Compared with con v entional A GC, it has faster con vergence speed. References [1] F . D ¨ orfler , J. W . Simpson-Porco, F . Bullo, Breaking the hierarchy: Dis- tributed control and economic optimality in microgrids, IEEE Trans. Con- trol Network Syst. 3 (3) (2016) 241–253. [2] Z. W ang, F . Liu, J. Z. Pang, S. H. Low , S. Mei, Distributed optimal fre- quency control considering a nonlinear network-preserving model, IEEE T rans. Power Syst. 34 (1) (2019) 76–86. [3] F . C. Schweppe, R. D. T abors, J. L. Kirtley , H. R. Outhred, F . H. Pickel, A. J. Cox, Homeostatic utility control, IEEE Trans. Power Apparatus Syst. P AS-99 (3) (1980) 1151–1163. [4] X. Zhang, A. Papachristodoulou, A real-time control framework for smart power networks: Design methodology and stability , Automatica 58 (2015) 43–50. 6 Bus1 Bus3 Bus4 Bus7 Bus8  0 20 40 60 80 100 120 140 160 180 Time (s) -1 -0.5 0 0.5 1 1.5 2 2.5 Bus19 Bus20 Bus21 Bus22 Bus23 0 20 40 60 80 100 120 140 160 180 Time (s) Figure 5: Dynamics of µ when line power congestion exists [5] T . Stegink, C. De Persis, A. van der Schaft, A unifying energy-based approach to stability of power grids with market dynamics, IEEE Trans. Autom. Control 62 (6) (2017) 2612–2622. [6] N. Li, C. Zhao, L. Chen, Connecting automatic generation control and economic dispatch from an optimization view , IEEE Trans. Control Net- work Syst. 3 (3) (2016) 254–264. [7] D. Cai, E. Mallada, A. W ierman, Distrib uted optimization decomposition for joint economic dispatch and frequenc y re gulation, IEEE Trans. Po wer Syst. 32 (6) (2017) 4370–4385. [8] C. Zhao, U. T opcu, N. Li, S. H.Low ., Design and stability of load-side primary frequency control in po wer systems, IEEE T rans. Autom. Control 59 (5) (2014) 1177–1189. [9] E. Mallada, C. Zhao, S. Low , Optimal load-side control for frequency regulation in smart grids, IEEE Trans. Autom. Control 62 (12) (2017) 6294–6309. [10] Z. W ang, F . Liu, S. H. Lo w, C. Zhao, S. Mei, Distributed frequency con- trol with operational constraints, part i: Per-node po wer balance, IEEE T rans. Smart Grid 10 (1) (2019) 40–52. [11] A. Kasis, E. Devane, C. Spanias, I. Lestas, Primary frequency regulation with load-side participation–part i: Stability and optimality , IEEE Trans. Power Syst. 32 (5) (2016) 3505–3518. [12] Z. W ang, F . Liu, S. H. Lo w, C. Zhao, S. Mei, Distributed frequency con- trol with operational constraints, part ii: Network power balance, IEEE T rans. Smart Grid 10 (1) (2019) 53–64. [13] C. Zhao, E. Mallada, S. H. Low , J. Bialek, Distributed plug-and-play op- timal generator and load control for power system frequency regulation, Int. J. Electr . Po wer Energy Syst. 101 (2018) 1–12. [14] Z. W ang, F . Liu, S. H. Low , P . Y ang, S. Mei, Distributed load-side control: Coping with variation of renewable generations, arXiv preprint arXiv:1804.04941 (2018). [15] X. Zeng, P . Yi, Y . Hong, L. Xie, Distributed continuous-time algorithms for nonsmooth extended monotropic optimization problems, SIAM J. Control Optim. 56 (6) (2018) 3973–3993. [16] F . H. Clarke, Optimization and nonsmooth analysis, V ol. 5, Siam, 1990. [17] H. Bauschke, P . L. Combettes, Conv ex Analysis and Monotone Operator Theory in Hilbert Spaces, Springer , 2017. [18] F . Facchinei, J.-S. Pang, Finite-dimensional variational inequalities and complementarity problems, Springer-V erlag, New Y ork, 2003. [19] Q. Liu, J. W ang, A one-layer projection neural network for nonsmooth op- timization subject to linear equalities and bound constraints, IEEE T rans. Neural Networks Learn. Syst. 24 (5) (2013) 812–824. [20] S. Boyd, L. V andenberghe, Conv ex optimization, Cambridge university press, 2004. [21] O. L. Mangasarian, R. Me yer , Nonlinear perturbation of linear programs, SIAM J. Control Optim. 17 (6) (1979) 745–752. [22] H. Zhou, X. Zeng, Y . Hong, Adapti v e exact penalty design for constrained distributed optimization, IEEE T rans. Autom. Control, in press (2019). [23] A. P . Ruszczy ´ nski, A. Ruszczynski, Nonlinear optimization, V ol. 13, Princeton univ ersity press, 2006. [24] J. Cortes, Discontinuous dynamical systems, IEEE Control Syst. Mag. 28 (3) (2008) 36–73. [25] K. W . Cheung, J. Cho w , G. Rogers, Po wer system toolbox ver . 3.0., Rens- selaer Polytechnic Institute and Cherry T ree Scientific Software (2009). 7