Distributed Optimal Load Frequency Control Considering Nonsmooth Cost Functions

# Distributed Optimal Load Frequency Control Considering Nonsmooth Cost Functions

Zhaojian Wang Feng Liu Changhong Zhao Zhiyuan Ma Wei Wei State Key Laboratory of Power Systems, Department of Electrical Engineering, Tsinghua University, Beijing 100084, China National Renewable Energy Laboratory, Golden, CO, 80401, US
###### Abstract

This work addresses the distributed frequency control problem in power systems considering controllable load with a nonsmooth cost. The nonsmoothness exists 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 flow are considered. Then, a distributed controller is derived using the Clark generalized gradient. We also prove 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 verifies the effectiveness of the proposed method.

###### keywords:
journal: xxx\setstretch

1

mytitlenotemytitlenotefootnotetext: This work was supported by the National Natural Science Foundation of China ( No. 51677100, U1766206, No. 51621065).

## 1 Introduction

With the proliferation of renewable generations, frequency control in power systems is facing a great challenge as power mismatch can fluctuate rapidly in a large amount. In this situation, the conventional centralized hierarchical control architecture may not respond fast enough due to large inertia of the traditional synchronous generators dorfler2016breaking (); wang2019distributed (). On the other hand, load-side controllable resources with fast response capabilities provide a new opportunity to frequency regulation Schweppe1980Homeostatic (). In addition, 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.

This work designs a distributed controller for the optimal load frequency control in power systems, where the regulation cost function can be nonsmooth. We relax the assumption of the objective function from being differentiable to nonsmooth. This work is partly motivated by zeng2018distributed (). However, different 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 difference is that the objective function in zeng2018distributed () is strictly convex with respect to all decision variables. That is not necessary in our work, where some variables may not appear in the objective function. In such a situation, we prove the asymptotic convergence 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. Section III formulates the optimal load frequency control problem and introduces the distributed controller. In Section IV, convergence of the closed-loop system and optimality of the equilibrium point are proved. We confirm the performance of the controller via simulations on IEEE 68-bus system in Section V. Section VI concludes the paper.

## 2 Problem Description

### 2.1 Preliminaries and notations

#### 2.1.1 Notations

In this paper, use () to denote the -dimensional (nonnegative) Euclidean space. For a column vector (matrix ), () denotes its transpose. For vectors , denotes the inner product of . denotes the Euclidean norm of . Use 1 to denote the vector with all elements. For a matrix , stands for the entry in the -th row and -th column of . Use to denote the Cartesian product of the sets . Given a collection of for in a certain set , denotes the column vector with a proper dimension, and as its components.

#### 2.1.2 Preliminaries

Let be a locally Lipschitz continuous function and denote its Clarke generalized gradient by (clarke:optimization, , Page 27). For a continuous strictly convex function , we have , where and .

Define the projection of onto a closed convex set as

 PΩ(x)=argminy∈Ω∥x−y∥ (1)

Use to denote the identity operator, i.e., , . Define . We have (bauschke2011convex, , Chapter 23.1).

A basic property of a projection is

 (x−PΩ(x))T(y−PΩ(x))≤0, ∀x∈Rn, y∈Ω (2)

Moreover, we also have (facchinei2003finite, , Theorem 1.5.5)

 (PΩ(x)−PΩ(y))T(x−y)≥∥PΩ(x)−PΩ(y)∥2 (3)

Define , and then is differentiable and convex with respect to (liu2013one, , Lemma 4). Moreover, we have

 V(x) =12∥PΩ(x)−PΩ(y)∥2 −(x−PΩ(x))T(PΩ(y)−PΩ(x)) (4) ≥12∥PΩ(x)−PΩ(y)∥2≥0 (5) ∇V(x) =PΩ(x)−PΩ(y) (6)

where the inequality is due to (2). From (2.1.2), holds only when .

### 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 , where is the set of buses and is the set of edges (transmission lines). Let denote the number of lines. The buses are divided into two types: generator buses, denoted by and load buses, denoted by . A generator bus contains a generator ( possibly with certain aggregate load). A load bus has only load with no generator. The graph is treated as directed with an arbitrary orientation and use or interchangeably to denote a directed edge from to . Without loss of generality, we assume the graph is connected and node is a reference node. The incidence matrix of the graph is denoted by , and we have .

We adopt a second-order linearized model to describe the frequency dynamics of each bus. We assume that the lines are lossless and adopt the DC power flow model Distributed_I:Wang (); Changhong:Design (). For each bus , let denote the rotor angle at node at time and the frequency. 111Sometimes, we also omit for simplicity. Let denote the controllable load. Let given constant denote any change in power injection, that occurs on the generation side or the load side, or both. Define as the angle difference between bus and , and its compact form is denoted by . Then for each node , the dynamics are

 ˙θij =ωi−ωj,j∈N (7a) ˙ωj =1Mj(Pmj−Plj−Djωj +∑i:i→jBijθij−∑k:j→kBjkθjk), j∈Ng (7b) 0 =Pmj−Plj−Djωj +∑i:i→jBijθij−∑k:j→kBjkθjk,j∈Nl (7c) where Mj>0 are inertia constants, Dj>0 are damping constants, and Bjk>0 are line parameters that depend on the reactance 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 . 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 propose a distributed controller based on the Clark generalized gradient to drive the power system to the optimal solution.

### 3.1 Optimization problem

The optimization problem is

 minPlj, ϕj f(Pl)=∑j∈Nfj(Plj) (8a) s.t. 0=Plj−Pmj−∑i:i→jBij(ϕi−ϕj) +∑k:j→kBjk(ϕj−ϕk), j∈N (8b) P––lj≤Plj≤¯¯¯¯Plj,j∈N (8c) θ–ij≤ϕi−ϕj≤¯¯¯θij,(i,j)∈E (8d)

where are constants, denoting the lower and upper bound of . are also constants, denoting the lower and upper bound of angle difference. The first constraint is the local power balance. is the virtual phase angle, which equals to at the optimal solution. Use to denote the virtual phase angle difference. In the DC power flow, we have , where is the power of line . Thus, (8d) is in fact the tie-line power limit constraint. We have the following assumptions.

###### Assumption 1.

is strictly convex.

###### Assumption 2.

The Slater’s condition (boyd2004convex, , Chapter 5.2.3) of (8) holds, i.e., problem (8) is feasible provided that the constraints are affine.

###### Remark 1.

Assumption 1 could be further relaxed, since a non-strictly convex function can be strictly convexified by using a nonlinear perturbation mangasarian1979nonlinear ().

###### Remark 2.

Problem (8) allows the cost function to be nonsmooth, which is required to be differentiable in the existing literature Changhong:Design (); mallada2017optimal (); Distributed_I:Wang (); Kasis:Primary1 (); Distributed_II:Wang (); zhao2018distributed (); wang2018distributed (). 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 is absent in the objective function of (8). It makes the paper not a trivial application of zeng2018distributed (), i.e., the objective function is not required to be strictly convex to all the decision variables. It makes the convergence proof more challenging.

###### Remark 3.

In the existing literature, the controller usually involves the projection of a gradient onto a convex set. If the objective function is nonsmooth, it becomes the projection of a subdifferential set onto a convex set. In this situation, the existence of trajectories is not guaranteed zeng2018distributed (); zhou2019adaptive (), which makes existing load frequency control methods inapplicable to the nonsmooth case.

### 3.2 Controller Design

To help the controller design, we make a modification on the problem (8).

 minPlj, ϕj f(Pl)=∑j∈Nfj(Plj)+12∑j∈Nz2j (9a) s.t. (9b)

where . For any feasible solution to (8), . Thus, (8) and (9) have same solutions.

Define the sets

 Ωj:={Plj | P––lj≤Plj≤¯¯¯¯Plj}, Ω=∏nj=1Ωj (10)

Then, we give the controller for each controllable load, which is denoted by OLC.

 ˙dj ∈{p:p=−dj+Plj+ωj−gj(Plj)−zj−μj, gj(Plj)∈∂fj(Plj)} (11a) ˙μj =Plj−Pmj−∑i:i→jBijϕij+∑k:j→kBjkϕjk (11b) ˙ϕj =∑i:i→jBij(μi−μj)−∑k:j→kBjk(μj−μk) −∑(i,j)∈Eη−ij+∑(j,k)∈Eη−jk+∑(i,j)∈Eη+ij−∑(j,k)∈Eη+jk +∑i:i→jBij(zi−zj)−∑k:j→kBjk(zj−zk) (11c) ˙φ+ij =−φ+ij+η+ij+ϕij−¯¯¯θij (11d) ˙φ−ij =−φ−ij+η−ij+θ–ij−ϕij (11e) Plj =PΩj(dj) (11f) η+ij =PR+(φ+ij) (11g) η−ij =PR+(φ−ij) (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 is difficult to measure. Similar to mallada2017optimal (); zhao2018distributed (); Distributed_II:Wang (), in (11b) can be substituted equivalently in following ways. For ,

 Plj−Pmj=−Mj˙ωj−Djωj+∑i:i→jPij−∑k:j→kPjk

For ,

 Plj−Pmj=−Djωj+∑i:i→jPij−∑k:j→kPjk

In this way, the measurement of load demand is avoided. We only need to measure , which are much easier to realize. Moreover, the power loss can be treated as unknown load demand, which can be also considered by this method.

## 4 Optimality and Convergence

In this section, we address the optimality of the equilibrium point and the convergence of the closed-loop system.

### 4.1 Optimality

Denote and . Let be an equilibrium of the closed-loop system (7), (11). Then, there exists such that

 0 =CTω∗ (12a) 0 =Pm−Pl∗−Dω∗−CBCTθ∗ (12b) 0 =−d∗+Pl∗+ω∗−g(Pl∗)−μ∗ (12c) 0 =Pl∗−Pm+CBCTϕ∗ (12d) 0 =−CBCTμ∗−Cη−∗+Cη+∗ (12e) 0 =−φ+∗+η+∗−CTϕ∗−¯¯¯θ (12f) 0 =−φ−∗+η−∗+θ–+CTϕ∗ (12g) Pl∗ =PΩ(d∗) (12h) η+∗ =PRm+(φ+∗) (12i) η−∗ =PRm+(φ−∗) (12j)

Now, we introduce the properties of the equilibrium points.

###### Theorem 1.

Suppose Assumptions 1 and 2 hold. We have

1. The nominal frequency is restored, i.e., for all .

2. If is an equilibrium point of (7), (11), then is an optimal solution to (8) and is an optimal solution to its dual problem.

3. for all . Moreover, the line limits are satisfied by , implying on every tie line .

4. At the equilibrium, is unique, with being unique up to (equilibrium) reference angles .

###### Proof.

1) From (12b) and (12d), we have . From (12a), we have with a constant . As is a diagonal positive definite matrix, we have .

2) From (12c) and (12f)-(12j), we have

 Pl∗ =PΩ(Pl∗−g(Pl∗)−μ∗) (13a) η+∗ =PRm+(η+∗−CTϕ∗−¯¯¯θ) (13b) η−∗ =PRm+(η−∗+θ–+CTϕ∗) (13c)

or equivalently,

 −g(Pl∗)−μ∗ ∈NΩ(Pl∗) (14a) −CTϕ∗−¯¯¯θ ∈NRm+(η+∗) (14b) θ–+CTϕ∗ ∈NRm+(η−∗) (14c)

By the KKT condition in (ruszczynski2006nonlinear, , 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 , which holds for any incidence matrix . Thus, we have with a constant . Then, we have . Moreover, by 2), we know , which implies that .

4) is unique because the objective function in (8a) is strictly convex in . is unique due to . By (12d), we know is unique modulo a rigid (uniform) rotation of all angles. Since , it implies that is also unique modulo a rigid rotation. This proves the uniqueness of . ∎

### 4.2 Convergence

Define the function

 V(x)=V1(x)+V2(x) (15)

where

 V1(x)=12∥∥Pl−Pl∗∥∥2+12∥μ−μ∗∥2+12∥ϕ−ϕ∗∥2 +12∥∥η+−η+∗∥∥2+12(θe−θ∗e)TB(θe−θ∗e) +12∥∥η−−η−∗∥∥2+12(ωg−ω∗g)TM(ωg−ω∗g) (16)
 V2(x)=−(d−Pl)T(Pl∗−Pl)−(φ+−η+)T(η+∗−η+) −(φ−−η−)T(η−∗−η−) (17)

Then, we have the following result about .

###### Lemma 2.

Suppose Assumptions 1 and 2 hold. Then the function has following properties

1. and holds only at the equilibrium point.

2. The time derivative of satisfies .

###### Proof.

1) By (2), we know that . From (15), (4.2) and (4.2), we know and holds only at the equilibrium point.

2) By (6), the gradient of is

 ∇V=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣∇dV∇μV∇ϕV∇η+V∇η−V∇θV∇ωgV⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣Pl−Pl∗μ−μ∗ϕ−ϕ∗η+−η+∗η−−η−∗B(θe−θ∗e)M(ωg−ω∗g)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (18)

Then, there is such that the time derivative of is

 ˙V=(Pl−Pl∗)T(−d+Pl+ω−g(Pl)−z−μ) +(μ−μ∗)T(Pl−Pm+CBCTϕ)+(θe−θ∗e)TBCTω +(ϕ−ϕ∗)T(−CBCTμ−Cη−+Cη+−CBCTz) +(η+−η+∗)T(−φ++η+−CTϕ−¯¯¯θ) +(η−−η−∗)T(−φ−+η−+θ–+CTϕ) +(ω−ω∗)T(Pm−Pl−Dω−CBθe) (19)

where the last item is due to the fact that, for each

 0 =(ωj−ω∗j)(Pmj−Plj−Djωj +∑i:i→jBijθij−∑k:j→kBjkθjk) (20)

Combing (4.2) and (12), we have

 ˙V =(~Pl)T(−~d+~Pl+~ω−g(Pl)+g(Pl∗)−~z−~μ) +~ϕT(−CBCT~μ−C~η−+C~η+−CBCT~z) +(~η−)T(−~φ−+~η−+CT~ϕ) +~θTBCT~ω+~ωT(−~Pl−D~ω−CB~θ) =−(Pl−Pl∗)T(d−d∗)+∥∥Pl−Pl∗∥∥2 (21a) (21b) (21c) −(Pl−Pl∗)T(g(Pl)−g(Pl∗))−~ωTD~ω (21d) −(~Pl)T~z−~ϕTCBCT~z (21e)

where .

By the (facchinei2003finite, , Theorem 1.5.5), we have

 −(Pl−Pl∗)T(d−d∗)+∥∥Pl−Pl∗∥∥2≤0 (22)

Similarly, and also hold.

The convexity of implies that

 −(Pl−Pl∗)T(g(Pl)−g(Pl∗))≤0 (23)

We also have because is positive definite. In addition,

 −(~Pl)T~z−~ϕTCBCT~z=−~zT⋅~z≤0 (24)

Then, (21a)-(21e) are all nonpositive, i.e., . ∎

The following result shows the stability of the closed-loop system (7), (11).

###### Theorem 3.

Suppose Assumptions 1 and 2 hold. Then the trajectory of the closed loop system (7), (11) has following properties

1. is bounded.

2. converges to equilibrium of the closed-loop system (7), (11).

3. The convergence of is a point, i.e., as for some equilibrium point .

###### Proof.

1) From Lemma 2, we know that is bounded. By (7c), is also bounded. Since is compact, there exists a constant such that

 ∥∥Pl+ω−g(Pl)−μ∥∥

Define following function

 ~Vd(d)=12∥d∥2 (26)

The time derivative of along the closed-loop system is

 ˙~Vd =dT(−d+Pl+ω−g(Pl)−μ) =−∥d∥2+dT(Pl+ω−g(Pl)−μ) ≤−∥d∥2+a1∥d∥=−2~Vd+a1√2~Vd (27)

Thus, is bounded, so is . Similarly, we can also have that and are bounded.

2) By the invariance principle in (cortes2008discontinuous, , Theorem 2), we know that the trajectory converges to the largest weakly invariant subset contained in , i.e., once a trajectory enters this subset, it will never departure from it.

From , we know , i.e., . Note that if due to the strict convexity of . Thus, we have in the set . Moreover, from (24), we have , or , which implies that . Then, we have from (11b). Similarly, by (11). Up to now, we know is constant except for .

Moreover, the equality in (3) holds only when or and . Thus, for , there are four combinations:

1. and ;

2. and , ;

3. and ;

4. and .

Thus, converges to equilibrium of the closed-loop system.

3) Fix any initial state and consider the trajectory of the closed-loop system. As is bounded, there exists an infinite sequence of time instants such that as , for some . Using this specific equlibrium point in the definition of , we have

 V∗=limt→∞V(x(t))=limtk→∞V(x(tk))=limx(tk)→^x∗V2(x(tk))=V2(^x∗)=0

Here, the first equality uses the fact that is nonincreasing in while lower-bounded, and therefore must render a limit value ; the second equality uses the fact that is the infinite subsequence of ; the third equality uses the fact that is absolutely continuous in ; the fourth equality is due to the continuity of , and the last equality holds as is an equilibrium point of .

The quadratic part implies that as . Moreover, from (12c), (12f) and (12g), we can get the corresponding . This completes the proof. ∎

## 5 Case studies

### 5.1 Test system

In this section, the IEEE 68-bus New England/New York interconnection test system Changhong:Design () is utilized to illustrate the performance of the proposed controller. The diagram of the 68-bus system is given in Fig.1. We run the simulation on Matlab using the Power System Toolbox PST (). Although the linear model is used in the analysis, the simulation model is much more detailed and realistic. The generator includes a two-axis subtransient reactance model, IEEE type DC1 exciter model, and a classical power system stabilizer model. AC (nonlinear) power flows are utilized, including non-zero line resistances. The upper bound of 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 objective function of each controllable load is

 f(Plj)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩(Plj)2−0.02,  Plj≤−0.212(Plj)2,  −0.2

It can be verified that is continuous, strictly convex and nonsmooth.

### 5.2 Simulation results

We consider the following scenario: at s, there is a step change of 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 power. In this subsection, we analyze the dynamic performance of the closed-loop system under the proposed controller OLC. In addition, automatic generation control (AGC) is tested in the same scenario as a benchmark. The setting of AGC is the same as that in zhao2018distributed (). The frequency dynamics under OLC and AGC are given in Fig.2. It is shown that both AGC and OLC can recover the frequency to the nominal value. They also have similar frequency nadir. Compared with AGC, the frequency under OLC has faster convergence speed.

The dynamics of and controllable load are illustrated in Fig.3. If tie-line power is not considered, then . From (12e), we know that , i.e., of each bus will converge to the same value. This is shown in the left of Fig.3, where of buses