Distributed Optimal Load Frequency Control Considering Nonsmooth Cost Functions

TL;DR

This paper develops a distributed control method for power system frequency regulation that effectively handles nonsmooth cost functions using the Clark generalized gradient, ensuring optimality and stability.

Contribution

It introduces a novel distributed control approach for nonsmooth cost functions in power systems, with proven optimality and stability.

Findings

Effective frequency regulation in IEEE 68-bus system

Handles nonsmooth cost functions with generalized gradient

Ensures asymptotic stability of the control system

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.

Figures

Equations

$${\mathcal{P}}_{\Omega }(x)=\arg {\min }_{y\in \Omega }\Vert x-y\Vert$$

$$(x-{\mathcal{P}}_{\Omega }(x){)}^{\mathrm{T}}(y-{\mathcal{P}}_{\Omega }(x))\le 0,\forall x\in {\mathbb{R}}^n,y\in \Omega$$

$$({\mathcal{P}}_{\Omega }(x)-{\mathcal{P}}_{\Omega }(y){)}^{\mathrm{T}}(x-y)\ge {\Vert {\mathcal{P}}_{\Omega }(x)-{\mathcal{P}}_{\Omega }(y)\Vert }^2$$

$$V(x)$$

$$-(x-{\mathcal{P}}_{\Omega }(x){)}^{\mathrm{T}}({\mathcal{P}}_{\Omega }(y)-{\mathcal{P}}_{\Omega }(x))$$

$$\ge \frac{1}{2}{\Vert {\mathcal{P}}_{\Omega }(x)-{\mathcal{P}}_{\Omega }(y)\Vert }^2\ge 0$$

$$\nabla V(x)$$

$${\dot{\theta }}_{ij}$$

$${\dot{\omega }}_j$$

$$\displaystyle\quad+\sum\nolimits_{i:i\rightarrow j}B_{ij}\theta_{ij}-\sum\nolimits_{k:j\rightarrow k}B_{jk}\theta_{jk}\bigg{)},\ j\in\mathcal{N}_{g}$$

$$0$$

$$+{\sum }_{i:i\to j}{B}_{ij}{\theta }_{ij}-{\sum }_{k:j\to k}{B}_{jk}{\theta }_{jk},j\in {\mathcal{N}}_l$$

$$\underset{P_j^l,{\varphi }_j}{\min }$$

$$0=P_j^l-P_j^m-{\sum }_{i:i\to j}{B}_{ij}({\varphi }_i-{\varphi }_j)$$

$$+{\sum }_{k:j\to k}{B}_{jk}({\varphi }_j-{\varphi }_k),j\in \mathcal{N}$$

$${\underline{P}}_j^l\le P_j^l\le {\overline{P}}_j^l,j\in \mathcal{N}$$

$${\underline{\theta }}_{ij}\le {\varphi }_i-{\varphi }_j\le {\overline{\theta }}_{ij},(i,j)\in \mathcal{E}$$

$$\underset{P_j^l,{\varphi }_j}{\min }$$

$$\text{eqref}OL{C}_2,\text{eqref}OL{C}_3,\text{eqref}OL{C}_4$$

$${\Omega }_j:=\left\{P_j^l|{\underline{P}}_j^l\le P_j^l\le {\overline{P}}_j^l\right\},\Omega ={\prod }_{j=1}^n{\Omega }_j$$

$${\dot{d}}_j$$

$$g_j(P_j^l)\in \partial f_j(P_j^l)\}$$

$${\dot{\mu }}_j$$

$${\dot{\varphi }}_j$$

$$-\sum _{(i,j)\in \mathcal{E}}{\eta }_{ij}^{-}+\sum _{(j,k)\in \mathcal{E}}{\eta }_{jk}^{-}+\sum _{(i,j)\in \mathcal{E}}{\eta }_{ij}^{+}-\sum _{(j,k)\in \mathcal{E}}{\eta }_{jk}^{+}$$

$$+\sum _{i:i\to j}{B}_{ij}(z_i-z_j)-\sum _{k:j\to k}{B}_{jk}(z_j-z_k)$$

$${\dot{\phi }}_{ij}^{+}$$

$${\dot{\phi }}_{ij}^{-}$$

$$P_j^l$$

$${\eta }_{ij}^{+}$$

$${\eta }_{ij}^{-}$$

$$P_j^l-P_j^m=-M_j{\dot{\omega }}_j-D_j{\omega }_j+{\sum }_{i:i\to j}{P}_{ij}-{\sum }_{k:j\to k}{P}_{jk}$$

$$P_j^l-P_j^m=-D_j{\omega }_j+{\sum }_{i:i\to j}{P}_{ij}-{\sum }_{k:j\to k}{P}_{jk}$$

$$0$$

$$P^{l\ast }$$

$${\eta }^{+\ast }$$

$${\eta }^{-\ast }$$

$$P^{l\ast }$$

$${\eta }^{+\ast }$$

$${\eta }^{-\ast }$$

$$-g(P^{l\ast })-{\mu }^{\ast }$$

Full text

††journal: xxxmytitlenotemytitlenotefootnotetext: 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 [1, 2]. On the other hand, load-side controllable resources with fast response capabilities provide a new opportunity to frequency regulation [3]. 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.

Recently, the so-called reverse engineering methodology is proposed by combining frequency control with optimal operation problems in power systems [4, 5, 6, 7, 8]. Under this framework, distributed load frequency control is widely investigated [8, 9, 10, 11, 12, 13, 14]. In [8], an optimal load frequency control problem is formulated and a distributed controller is derived using controllable loads to realize primary frequency control. To eliminated the frequency deviation, the method is further extended in [9, 13] to realize a secondary load frequency control. At the same time, the tie-line power limit is considered. The design approach is generalized in [11], where the specific model requirement is eliminated. It only requires that the bus dynamics satisfy a passivity condition to guarantee asymptotic stability. In [10, 12], the operational constraints including regulation capacity limits and tie-line power limits are considered, which guarantee both steady-state and transient capacity limit constraints. In [14], the distributed load frequency control under time-varying and unknown power injection is investigated, which can recover the nominal frequency even 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 power. In the existing literature, the cost of controllable load is assumed to be differentiable, or equivalently, 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 different intervals. In such a situation, the regulation is inherently nonsmooth, which makes existing methods difficult to apply.

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 [15]. 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 [15] 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 $\mathbb{R}^{n}$ ($\mathbb{R}^{n}_{+}$) to denote the $n$-dimensional (nonnegative) Euclidean space. For a column vector $x\in\mathbb{R}^{n}$ (matrix $A\in\mathbb{R}^{m\times n}$), $x^{\mathrm{T}}$($A^{\mathrm{T}}$) denotes its transpose. For vectors $x,y\in\mathbb{R}^{n}$, $x^{\mathrm{T}}y=\left\langle x,y\right\rangle$ denotes the inner product of $x,y$. $\left\|x\right\|=\sqrt{x^{\mathrm{T}}x}$ denotes the Euclidean norm of $x$. Use 1 to denote the vector with all $1$ elements. For a matrix $A=[a_{ij}]$, $a_{ij}$ stands for the entry in the $i$-th row and $j$-th column of $A$. Use $\prod_{i=1}^{n}\Omega_{i}$ to denote the Cartesian product of the sets $\Omega_{i},i=1,\cdots,n$. Given a collection of $y_{i}$ for $i$ in a certain set $Y$, $y$ denotes the column vector $y:=(y_{i},i\in Y)$ with a proper dimension, and $y_{i}$ as its components.

2.1.2 Preliminaries

Let $f(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$ be a locally Lipschitz continuous function and denote its Clarke generalized gradient by $\partial f(x)$ [16, Page 27]. For a continuous strictly convex function $f(x):\mathbb{R}^{n}\rightarrow\mathbb{R}$, we have $(g_{x}-g_{y})^{\mathrm{T}}(x-y)>0,\ \forall x\neq y$, where $g_{x}\in\partial f(x)$ and $g_{y}\in\partial f(y)$.

Define the projection of $x$ onto a closed convex set $\Omega$ as

[TABLE]

Use ${\rm{Id}}$ to denote the identity operator, i.e., ${\rm{Id}}(x)=x$, $\forall x$. Define $N_{\Omega}(x)=\{v|\left\langle v,y-x\right\rangle\leq 0,\forall y\in\Omega\}$. We have $\mathcal{P}_{\Omega}(x)=({\rm{Id}}+N_{\Omega})^{-1}(x)$ [17, Chapter 23.1].

A basic property of a projection is

[TABLE]

Moreover, we also have [18, Theorem 1.5.5]

[TABLE]

Define V(x):=\frac{1}{2}\big{(}\left\|x-\mathcal{P}_{\Omega}(y)\right\|^{2}-\left\|x-\mathcal{P}_{\Omega}(x)\right\|^{2}\big{)}, and then $V(x)$ is differentiable and convex with respect to $x$ [19, Lemma 4]. Moreover, we have

[TABLE]

where the inequality is due to (2). From (2.1.2), $V(x)=0$ holds only when $\mathcal{P}_{\Omega}(x)=\mathcal{P}_{\Omega}(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 $\mathcal{G}:=(\mathcal{N},\mathcal{E})$, where $\mathcal{N}=\{0,1,2,...n\}$ is the set of buses and $\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ is the set of edges (transmission lines). Let $m=|\mathcal{E}|$ denote the number of lines. The buses are divided into two types: generator buses, denoted by $\mathcal{N}_{g}$ and load buses, denoted by $\mathcal{N}_{l}$. A generator bus contains a generator ( possibly with certain aggregate load). A load bus has only load with no generator. The graph $\mathcal{G}$ is treated as directed with an arbitrary orientation and use $(i,j)\in\mathcal{E}$ or $i\rightarrow j$ interchangeably to denote a directed edge from $i$ to $j$. Without loss of generality, we assume the graph is connected and node [math] is a reference node. The incidence matrix of the graph is denoted by $C$, and we have $\textbf{1}^{\rm T}C=0$.

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 [10, 8]. For each bus $j\in\mathcal{N}$, let $\theta_{j}(t)$ denote the rotor angle at node $j$ at time $t$ and $\omega_{j}(t)$ the frequency. 111Sometimes, we also omit $t$ for simplicity. Let $P^{l}_{j}(t)$ denote the controllable load. Let given constant $P^{m}_{j}$ denote any change in power injection, that occurs on the generation side or the load side, or both. Define $\theta_{ij}=\theta_{i}-\theta_{j}$ as the angle difference between bus $i$ and $j$, and its compact form is denoted by $\theta_{e}=(\theta_{ij},(i,j)\in\mathcal{E})$. Then for each node $j\in\mathcal{N}$, the dynamics are

[TABLE]

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_{j}^{m}$. 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

[TABLE]

where $\underline{P}_{j}^{l}\leq\overline{P}_{j}^{l}$ are constants, denoting the lower and upper bound of $P_{j}^{l}$. ${\underline{\theta}}_{ij}\leq{\overline{\theta}}_{ij}$ are also constants, denoting the lower and upper bound of angle difference. The first constraint is the local power balance. $\phi_{j}$ is the virtual phase angle, which equals to $\theta_{j}$ at the optimal solution. Use $\phi_{ij}=\phi_{i}-\phi_{j}$ to denote the virtual phase angle difference. In the DC power flow, we have $P_{ij}=B_{ij}\theta_{ij}$, where $P_{ij}$ is the power of line $(i,j)$. Thus, (8d) is in fact the tie-line power limit constraint. We have the following assumptions.

Assumption** 1****.**

$f_{j}(P_{j}^{l})$ is strictly convex.

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 affine.

Remark 1**.**

Assumption 1 could be further relaxed, since a non-strictly convex function can be strictly convexified by using a nonlinear perturbation [21].

Remark 2**.**

Problem (8) allows the cost function $f_{j}(P_{j}^{l})$ to be nonsmooth, which is required to be differentiable 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 $\phi_{j}$ is absent in the objective function of (8). It makes the paper not a trivial application of [15], 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 [15, 22], 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).

[TABLE]

where $z_{j}=P_{j}^{l}-P_{j}^{m}-\sum\nolimits_{i:i\rightarrow j}B_{ij}\phi_{ij}+\sum\nolimits_{k:j\rightarrow k}B_{jk}\phi_{jk}$. For any feasible solution to (8), $z_{j}=0$. Thus, (8) and (9) have same solutions.

Define the sets

[TABLE]

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

[TABLE]

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_{j}^{l}$ is difficult to measure. Similar to [9, 13, 12], $P_{j}^{l}-P_{j}^{m}$ in (11b) can be substituted equivalently in following ways. For $j\in\mathcal{N}_{g}$,

[TABLE]

For $j\in\mathcal{N}_{l}$,

[TABLE]

In this way, the measurement of load demand $P_{j}^{l}$ is avoided. We only need to measure $\omega_{j},P_{ij}$, 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 $x=(\theta,\omega_{g},d,\mu,\phi,\varphi^{-},\varphi^{+})$ and $y=(x,P^{l},{\eta}^{+},{\eta}^{-})$. Let $x^{*}=(\theta^{*},\omega_{g}^{*},d^{*},\mu^{*},$ $\phi^{*},$ $\varphi^{-*},\varphi^{+*})$ be an equilibrium of the closed-loop system (7), (11). Then, there exists $g(P^{l*})\in\partial f(P^{l*})$ such that

[TABLE]

Now, we introduce the properties of the equilibrium points.

Theorem 1**.**

Suppose Assumptions 1 and 2 hold. We have

The nominal frequency is restored, i.e., $\omega^{*}_{j}=0$ for all $j\in\mathcal{N}$. 2. 2.

If $x^{*}$ is an equilibrium point of (7), (11), then $(P^{l*},\phi^{*})$ is an optimal solution to (8) and $(\mu^{*},{\eta}^{+*},{\eta}^{-*})$ is an optimal solution to its dual problem. 3. 3.

$\phi^{*}_{ij}=\theta^{*}_{ij}$ for all $(i,j)\in\mathcal{E}$. Moreover, the line limits are satisfied by $x^{*}$, implying $\underline{\theta}_{ij}\leq\theta^{*}_{ij}\leq\overline{\theta}_{ij}$ on every tie line $(i,j)\in\mathcal{E}$. 4. 4.

At the equilibrium, $(\theta^{*},\phi^{*},\omega_{g}^{*},P^{l*})$ is unique, with $(\theta^{*},\phi^{*})$ being unique up to (equilibrium) reference angles $(\theta_{0},\phi_{0})$.

Proof.

1. From (12b) and (12d), we have $\textbf{1}^{\rm T}D\omega^{*}=0$. From (12a), we have $\omega^{*}=\omega_{0}\cdot\textbf{1}$ with a constant $\omega_{0}$. As $D$ is a diagonal positive definite matrix, we have $\omega_{0}=0$.

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

[TABLE]

or equivalently,

[TABLE]

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 $CBC^{\rm T}(\theta^{*}-\phi^{*})=0$, which holds for any incidence matrix $C$. Thus, we have $\theta^{*}-\phi^{*}=c_{0}\cdot\textbf{1}$ with a constant $c_{0}$. Then, we have $\theta_{ij}^{*}-\phi_{ij}^{*}=0$. Moreover, by 2), we know $\underline{\theta}_{ij}\leq\phi_{ij}^{*}\leq\overline{\theta}_{ij}$, which implies that $\underline{\theta}_{ij}\leq\theta_{ij}^{*}\leq\overline{\theta}_{ij}$.

4. $P^{l*}$ is unique because the objective function in (8a) is strictly convex in $P^{l}$. $\omega^{*}$ is unique due to $\omega^{*}=0$. By (12d), we know $\phi^{*}$ is unique modulo a rigid (uniform) rotation of all angles. Since $\theta^{*}-\phi^{*}=c_{0}\cdot\textbf{1}$, it implies that $\theta^{*}$ is also unique modulo a rigid rotation. This proves the uniqueness of $(\theta^{*},\phi^{*},\omega_{g}^{*},P^{l*})$. ∎

4.2 Convergence

Define the function

[TABLE]

where

[TABLE]

[TABLE]

Then, we have the following result about $V(x)$.

Lemma 2**.**

Suppose Assumptions 1 and 2 hold. Then the function $V(x)$ has following properties

$V(x)\geq 0$ and $V(x)=0$ holds only at the equilibrium point. 2. 2.

The time derivative of $V(x(t))$ satisfies $\dot{V}(x(t))\leq 0$.

Proof.

1. By (2), we know that $V_{2}(x)\geq 0$. From (15), (4.2) and (4.2), we know $V(x)\geq 0$ and $V(x)=0$ holds only at the equilibrium point.

2. By (6), the gradient of $V$ is

[TABLE]

Then, there is $g(P^{l})\in\partial f(P^{l})$ such that the time derivative of $V$ is

[TABLE]

where the last item is due to the fact that, for each $j\in\mathcal{N}_{l}$

[TABLE]

Combing (4.2) and (12), we have

[TABLE]

where $\tilde{x}=x-x^{*}$.

By the [18, Theorem 1.5.5], we have

[TABLE]

Similarly, $-(\eta^{+}-\eta^{+*})^{\rm T}({\varphi}^{+}-{\varphi}^{+*})+\left\|\eta^{+}-\eta^{+*}\right\|^{2}\leq 0$ and $-(\eta^{-}-\eta^{-*})^{\rm T}({\varphi}^{-}-{\varphi}^{-*})+\left\|\eta^{-}-\eta^{-*}\right\|^{2}\leq 0$ also hold.

The convexity of $f$ implies that

[TABLE]

We also have $-\tilde{\omega}^{\rm T}D\tilde{\omega}\leq 0$ because $D$ is positive definite. In addition,

[TABLE]

Then, (21a)-(21e) are all nonpositive, i.e., $\dot{V}(x(t))\leq 0$. ∎

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

$(x(t),P^{l}(t),\eta^{+}(t),\eta^{-}(t))$ is bounded. 2. 2.

$(x(t),P^{l}(t),\eta^{+}(t),\eta^{-}(t))$ converges to equilibrium of the closed-loop system (7), (11). 3. 3.

The convergence of $x(t)$ is a point, i.e., $x(t)\rightarrow x^{*}$ as $t\rightarrow\infty$ for some equilibrium point $x^{*}$.

Proof.

From Lemma 2, we know that $(\theta(t),\omega_{g}(t),\mu(t),$ $\phi(t),$ $P^{l}(t),\eta^{+}(t),\eta^{-}(t))$ is bounded. By (7c), $\omega_{l}(t)$ is also bounded. Since $\partial f(P^{l})$ is compact, there exists a constant $a_{1}$ such that

[TABLE]

Define following function

[TABLE]

The time derivative of $\tilde{V}_{d}(d)$ along the closed-loop system is

[TABLE]

Thus, $\tilde{V}_{d}(d(t)),t\geq 0$ is bounded, so is $d(t),t\geq 0$. Similarly, we can also have that $\varphi^{+}(t),t\geq 0$ and $\varphi^{-}(t),t\geq 0$ are bounded.

2. By the invariance principle in [24, Theorem 2], we know that the trajectory $x(t)$ converges to the largest weakly invariant subset $W^{*}$ contained in $W:=\{\ x\ |\ \dot{V}(x)=0\ \}$, i.e., once a trajectory enters this subset, it will never departure from it.

From $\tilde{\omega}^{\rm T}D\tilde{\omega}=0$, we know $\tilde{\omega}=0$, i.e., $\omega(t)=\omega^{*}$. Note that $(P^{l}-P^{l*})^{\rm T}\left(g(P^{l})-g(P^{l*})\right)>0$ if $x\neq x^{*}$ due to the strict convexity of $f(P^{l})$. Thus, we have $P^{l}(t)=P^{l*}$ in the set $W^{*}$. Moreover, from (24), we have $\tilde{P}^{l}(t)=CBC^{\rm T}\tilde{\phi}(t)$, or $0=\dot{\tilde{P}}^{l}(t)=CBC^{\rm T}\dot{\tilde{\phi}}(t)$, which implies that $\dot{\tilde{\phi}}(t)=0$. Then, we have $\dot{\mu}_{j}(t)=0$ from (11b). Similarly, $\dot{d}_{j}(t)=0$ by (11). Up to now, we know $x(t)$ is constant except $\varphi^{-}(t),\varphi^{+}(t)$ for $t\rightarrow\infty$.

Moreover, the equality in (3) holds only when $\mathcal{P}_{\Omega}(x)=\mathcal{P}_{\Omega}(y)$ or $x=\mathcal{P}_{\Omega}(x)$ and $y=\mathcal{P}_{\Omega}(y)$. Thus, for $\eta^{+}(t),\ \eta^{-}(t),t\rightarrow\infty$, there are four combinations:

$\eta^{+}(t)=\eta^{+*}$ and $\eta^{-}(t)=\eta^{-*}$; 2. 2.

$\eta^{+}(t)=\eta^{+*}$ and ${\varphi}^{-}(t)=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{-}(t)\right)=\eta^{-}(t)$, ${\varphi}^{-*}=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{-*}\right)=\eta^{-*}$; 3. 3.

${\varphi}^{+}(t)=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{+}(t)\right)=\eta^{+}(t),{\varphi}^{+*}=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{+*}\right)=\eta^{+*}$ and $\eta^{-}=\eta^{-*}$; 4. 4.

${\varphi}^{+}(t)=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{+}(t)\right)=\eta^{+}(t),{\varphi}^{+*}=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{+*}\right)=\eta^{+*}$ and ${\varphi}^{-}(t)=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{-}(t)\right)=\eta^{-}(t),{\varphi}^{-*}=\mathcal{P}_{\mathbb{R}^{m}_{+}}\left({\varphi}^{-*}\right)=\eta^{-*}$.

Thus, $(x(t),P^{l}(t),\eta^{+}(t),\eta^{-}(t))$ converges to equilibrium of the closed-loop system.

3. Fix any initial state $x(0)$ and consider the trajectory $(x(t),t\geq 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})\to\hat{x}^{*}$ as $t_{k}\to\infty$, for some $\hat{x}^{*}\in W^{*}$. Using this specific equlibrium point $\hat{x}^{*}$ in the definition of $V$, we have

[TABLE]

Here, the first equality uses the fact that $V(t)$ is nonincreasing in $t$ while lower-bounded, and therefore must render a limit value $V^{*}$; the second equality uses the fact that $t_{k}$ is the infinite subsequence 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 $\hat{x}^{*}$ is an equilibrium point of $V$.

The quadratic part $V_{1}$ implies that $(\theta,\omega_{g},P^{l},\mu,\phi,\eta^{-},$ $\eta^{+})\to(\theta^{*},\omega_{g}^{*},P^{l*},\mu^{*},\phi^{*},\eta^{-*},\eta^{+*})$ as $t\to\infty$. Moreover, from (12c), (12f) and (12g), we can get the corresponding $d^{*},\varphi^{-*},\varphi^{+*}$. 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 [8] 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 [25]. 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 $P_{j}^{l}$ 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

[TABLE]

It can be verified that $f(P_{j}^{l})$ is continuous, strictly convex and nonsmooth.

5.2 Simulation results

We consider the following scenario: at $t=1$s, 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 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 [13]. 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 $\mu$ and controllable load are illustrated in Fig.3. If tie-line power is not considered, then $\eta^{-*}=\eta^{+*}=0$. From (12e), we know that $CBC^{\rm T}\mu^{*}=0$, i.e., $\mu$ of each bus will converge to the same value. This is shown in the left of Fig.3, where $\mu$ of buses $1\sim 5$ is given. The right part of Fig.3 illustrates the dynamics of controllable load at these buses. They also converge to the same value except bus 2 and bus 5 as the objective 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.

We further take the tie-line limits into account. In this scenario, the active power limit of line $(1,2)$ is set to be [math]. Then, its tie-line power dynamics with and without limit is given in Fig.4. The active power is $-2.1$p.u. if there is no limit. It decreases to zero if the limit is considered, which verifies the effectiveness of the line power control.

If some tie-line limit is reached, the $\mu$ of buses near from this line will diverge. This is illustrated in the left part of Fig.5, where $\mu_{1},\mu_{3},\mu_{4}$ are all different. We also find that the $\mu$ of buses far from the line still converges to the same value, which is given in the right part of Fig.5. It is shown that $\mu_{19}\sim\mu_{23}$ converge to $0.83$.

6 Conclusion

In this paper, we investigate the distributed load frequency control in power systems when the regulation cost is nonsmooth. 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 address 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 proved that the convergence is to a specific point. Finally, numerical experiments on the IEEE 68-bus system show that the frequency is recovered to the nominal value. Compared with conventional AGC, it has faster convergence speed.