Robust Distributed Load Frequency Control for Multi-Area Power Systems with Photovoltaic and Battery Energy Storage System

Abstract

The intermittent power generation of renewable energy sources (RESs) interrupts the balance between power generation and demand load due to the increased frequency fluctuation, which challenges the frequency stability analysis and control synthesis of power generation systems. This paper proposes a robust distributed load frequency control (DLFC) scheme for multi-area power systems. Firstly, a multi-area power system is constructed by integrating photovoltaic (PV) and battery energy storage systems (BESSs). Then, by employing the linear matrix inequality (LMI) technique, the sufficient condition capable of ensuring that the proposed controller satisfies $H_{\infty }$ robust performance in the sense of asymptotic stability is derived. Finally, testing is conducted on a four-area renewable power system, and results verify the strong robustness of the proposed controller against load disturbance and intermittence of RESs.

1. Introduction

Renewable energy sources (RESs) have been experiencing significant growth in recent decades, due to sustainable environmental development and decarbonization. It is estimated that the year U.S.’s 2050 end-use all-purpose load will be met with 42.5% photovoltaic (PV) and concentrated solar power with storage [1]. A battery energy storage system (BESS) is used to store electric energy in a battery. With the development of sealed, recombinant lead–acid battery technology and advanced compounds, the potential of large-scale grid applications of the BESS has been extended in recent decades [2]. A BESS is always coupled with RESs and used as a buffer to enable maximum output power and to have a smoother output profile [3]. By charging and discharging as a local energy storage system when frequency deviation appears, a BESS is a means of mitigating the intermittence of RESs [4,5]. Thus, the grid integration of RESs into power systems is reasonable, following the development of RESs. Though the control techniques for traditional fossil fuel power generation to maintain frequency stability have been well developed, the intermittent nature of RESs, the stochasticity of load disturbance, and the continuously changing power demand causes a power imbalance between generation and load demand [6], causing a deviation in the system frequency. Therefore, load frequency control (LFC) has garnered considerable attention in power systems over the past few decades.

The goal of LFC is to maintain the frequency at the expected value by regulating the frequency reference set point. Distributed load frequency control (DLFC) is widely employed in multi-area power systems, considering the interactions between the local control area and its neighbors. Though a centralized control strategy has better control performance, expanded computational and capacity complexities will lead to a heavy communication burden, due to frequent data exchange between areas [7]. The impeded communication may deteriorate the dynamic performance of control, especially when the number of involved areas is large [8]. Compared with centralized LFC, DLFC can save communication resources and have faster dynamic performance [9]. Due to its strengths of balancing communication resource allocation effectively and adjusting the frequency accurately, DLFC is more commonly adopted for multi-area power systems, rather than centralized and decentralized control strategies [9,10].

Proportional integral (PI) and proportional integral derivative (PID) controllers have been employed in LFC typically and traditionally. But control gains of these controllers are fixed and designed at nominal operating points, without considering the adaptation of system uncertainty and load disturbance [11]. Techniques like the internal model control (IMC) method [12], Particle Swarm Optimization, and the Genetic Algorithm [13] were proposed for tuning the parameters of the PI/PID controller. However, the tuning process of IMC degrades the closed-loop response and robustness [12]. And for PSO and GA, the weight function for optimization is hard to design with the stability guarantees. Other control techniques like sliding mode control [14], model predictive control [15,16], robust control [17], and fuzzy control [18] have also been proposed for LFC synthesis. Though sliding mode control is robust against matched load disturbances and parameter uncertainties, chattering is an inherent problem which can increase vibration in the frequency regulation [8]. Model predictive control has excellent tracking ability and capability to cope with uncertainties [4,19], but the terminal constraints are difficult to handle for ensuring stability [20]. Fuzzy control does well in coping with parameter uncertainty in system modeling but has limitations in inferring the reasonable fuzzy logic [8]. Among the aforementioned discussions, robust performance should be considered in the control synthesis to enhance the robustness of the systems with uncertainties [19]. In [21] and [22], the authors demonstrated that the robust LFC method has good capability against load disturbance and is able to deal with time-delay attack problems [23,24]. Though load disturbance and time-delay attack problems are important when designing the controller, the integration of RESs is of urgent need in future power systems from the sustainable development perspective. The works in [16,18,19] did not consider RESs in power systems. On the other hand, the case study in [24] only considered a single-area power system, and their control strategy had limitations in practical multi-area power systems. Though hydro power generation was employed in [25], hydro power systems are used as the main generation in one area; thus, there is still potential for exploring other categories of RESs. To sum up, some of the above studies were limited to the absence of RESs, while some were limited to one area; achieving a robust performance for LFC should be in consideration.

To address the limitations in the above literature, this work proposes a robust DLFC scheme for multi-area power systems with PV and BESSs. It can mitigate the frequency fluctuation caused by load disturbances and the intermittence of RESs. The main contributions of this work are summarized as follows:

• A novel distributed load frequency control scheme is proposed for a renewable power system comprising PV and BESSs. The proposed controller is designed to mitigate the impact of load disturbance and the intermittence of PV power.
• Robustness is guaranteed for the proposed load frequency control scheme by incorporating the $H_{\infty }$ robust performance index $\gamma$. Two theorems on sufficient conditions for linear matrix inequalities (LMIs) are derived to meet this objective.

The rest of this paper is organized as follows. Section 2 presents the dynamic models of each part of the power system. The proposed DLFC controller satisfying $H_{\infty }$ robust performance is described in Section 3. Section 4 presents case studies validating the effectiveness of the proposed robust controller, compared with the traditional PID controller. Section 5 presents sensitivity analysis of the impact of the capacity of PV power and the impact of load disturbance. Section 6 concludes the paper.

2. Modeling of the Renewable Power System

In this paper, the LFC model of multi-area interconnected power systems including PV and BESSs is constructed; the framework is shown in Figure 1. Index $i\in \left\{1, 2, \dots ,N\right\}$ indicates the ith area.

2.1. Modeling of PV

The PV output power (W) is expressed as [2,26]

$$P_{PV}=\eta S\mathsf{\Phi }\left\{1-0.005\left(T_a+25\right)\right\}$$

(1)

where $\eta$ is the conversion efficiency of the PV array and the range is 9% to 12%, $S$ is the measured area of the PV array (${\mathrm{m}}^2$), $\mathsf{\Phi }$ is the radiation intensity ($\mathrm{W}/{\mathrm{m}}^2$), and $T_a$ is ambient temperature (°C). Basically, $\eta$ and $S$ are constant, and the value of $P_{PV}$ depends on $T_a$ and $\mathsf{\Phi }$. In this paper, we assume that the ambient temperature is 25 °C. Thus, the PV output power is linearly proportional to $\mathsf{\Phi }$.

The transfer function of PV can be expressed as

$$G_{{pv}_i}\left(s\right)={\displaystyle \frac{K_{{PV}_i}}{1+s{T}_{{PV}_i}}}={\displaystyle \frac{\mathsf{\Delta }{P}_{{PV}_i}}{\mathsf{\Delta }{\mathsf{\Phi }}_i}}$$

(2)

where $K_{{PV}_i}$ is the PV gain factor, $T_{{PV}_i}$ is the PV time constant, and $∆P_{PV}$ is the PV output power deviation.

2.2. Modeling of BESS

An energy storage system is used to supply insufficient energy in this renewable power system to help maintain system stability. The transfer function of the BESS can be expressed as [2]

$$G_{{BESS}_i}\left(S\right)={\displaystyle \frac{K_{{BESS}_i}}{1+s{T}_{{BESS}_i}}}={\displaystyle \frac{\mathsf{\Delta }{P}_{{BESS}_i}}{\mathsf{\Delta }{f}_i}}$$

(3)

The dynamic of the BESS can be expressed as

$${\displaystyle \frac{d\mathsf{\Delta }{P}_{{BESS}_i}}{dt}}={\displaystyle \frac{K_{{BESS}_i}}{T_{{BESS}_i}}}\mathsf{\Delta }{f}_i-{\displaystyle \frac{1}{T_{{BESS}_i}}}\mathsf{\Delta }{P}_{{BESS}_i}$$

(4)

where $K_{{BESS}_i}$ is the BESS gain factor, $T_{{BESS}_i}$ is the BESS time constant, $\mathsf{\Delta }{P}_{{BESS}_i}$ is the deviation of power provided by BESS, and $\mathsf{\Delta }{f}_i$ is the frequency deviation of area i.

2.3. Modeling of LFC

According to Figure 1, the dynamics of the generator can be described as

$${\displaystyle \frac{d\mathsf{\Delta }{P}_{m_i}}{dt}}=-{\displaystyle \frac{1}{T_{chi}}}\mathsf{\Delta }{P}_{m_i}+{\displaystyle \frac{1}{T_{chi}}}\mathsf{\Delta }{P}_{v_i}$$

(5)

where $T_{chi}$ is the time constant of the generator, $\mathsf{\Delta }{P}_{v_i}$ is the generator value position deviation, and $\mathsf{\Delta }{P}_{m_i}$ is the generator mechanical power deviation.

The generator value position deviation is determined by the frequency deviation and control signal. Its dynamics can be described as

$${\displaystyle \frac{d\mathsf{\Delta }{P}_{v_i}}{dt}}=-{\displaystyle \frac{1}{T_{g_i}}}\mathsf{\Delta }{P}_{v_i}-{\displaystyle \frac{1}{R_i{T}_{g_i}}}\mathsf{\Delta }{f}_i+{\displaystyle \frac{1}{T_{g_i}}}{\mathsf{\Delta }P}_{c_i}$$

(6)

where $T_{g_i}$ is the time constant of the governor, ${\mathsf{\Delta }P}_{c_i}$ is the deviation of load reference set-point and $R_i$ is the speed drop coefficient.

Ignoring the power loss in transmission of the tie-line, the dynamics of tie-line active power deviation can be expressed as

$${\displaystyle \frac{d\mathsf{\Delta }{P}_{tie}^i}{dt}}=2\pi \sum _{j=1, j\ne i}^N{T}_{ij}\left(\mathsf{\Delta }{f}_i-\mathsf{\Delta }{f}_j\right)$$

(7)

where $T_{ij}$ is the tie-line synchronizing coefficient between area i and area j, which is used to quantify the coupling between two interconnected areas (area i and area j) in terms of power exchange through the tie-line.

LFC is to control the frequency deviation and tie-line power fluctuation within the limited values; the output should be

$${ACE}_i={\beta }_i\mathsf{\Delta }{f}_i+\mathsf{\Delta }{P}_{tie}^i$$

(8)

where ${ACE}_i$ is the area control error, ${\beta }_i$ is the frequency bias constant, ${\beta }_i={\displaystyle \frac{1}{R_i}}+D_i$, $\mathsf{\Delta }{P}_{tie}^i$ is the tie-line active power deviation between area i and other neighbor areas.

The dynamic relationship between the incremental mismatch power $(∆P_{m_i}+P_{{BESS}_i}-∆P_{tie}^i-∆P_{PV}-∆P_{L_i})$ and the frequency deviation $(∆f_i)$ can be described as

$${\displaystyle \frac{d\mathsf{\Delta }{f}_i}{dt}}=-{\displaystyle \frac{D_i}{T_{m_i}}}\mathsf{\Delta }{f}_i+{\displaystyle \frac{1}{T_{m_i}}}\left(\mathsf{\Delta }{P}_{m_i}+\mathsf{\Delta }{P}_{{BESS}_i}-\mathsf{\Delta }{P}_{tie}^i-∆P_{{PV}_i}-∆P_{L_i}\right)$$

(9)

where $D_i$ is the equivalent damping coefficient of the generator, ${T_m}_i$ is the equivalent inertia of the generator, $∆P_{L_i}$ is the load fluctuation.

Based on (5)–(9), the system state can be expressed as $x_i={\left[{\mathsf{\Delta }f}_i {\mathsf{\Delta }P}_{m_i} \mathsf{\Delta }{P}_{v_i} \mathsf{\Delta }{P}_{tie}^i \mathsf{\Delta }{P}_{{BESS}_i} \int {ACE}_i\right]}^T$. The system dynamics of area i is depicted as

$$\left\{\begin{array}{l}{\dot{x}}_i\left(t\right)={\mathcal{A}}_i{x}_i\left(t\right)+{\mathcal{B}}_i{u}_i\left(t\right)+{\displaystyle \sum _{j=1,j\ne i}^N}{\mathcal{A}}_{ij}{x}_j\left(t\right)+{\mathcal{F}}_i\left(\mathsf{\Delta }{P}_{{PV}_i}\left(t\right)+\mathsf{\Delta }{P}_{L_i}\left(t\right)\right)\\ y_i\left(t\right)={\mathcal{C}}_i{x}_i\left(t\right) \\ z_i\left(t\right)={\mathcal{D}}_i{x}_i\left(t\right) \end{array}\right.$$

(10)

where $u_i$ is the control input, $y_i$ is the measured output, and $z_i$ indicates the interested signal,

$$\begin{array}{lll}\hfill {\mathcal{A}}_i& =& \left[\begin{array}{cccccc}-{\displaystyle \frac{D_i}{T_{m_i}}}& {\displaystyle \frac{1}{T_{m_i}}}& 0& -{\displaystyle \frac{1}{T_{m_i}}}& {\displaystyle \frac{1}{T_{m_i}}}& 0\\ 0& -{\displaystyle \frac{1}{T_{chi}}}& {\displaystyle \frac{1}{T_{chi}}}& 0& 0& 0\\ -{\displaystyle \frac{1}{R_i{T}_{g_i}}}& 0& -{\displaystyle \frac{1}{T_{g_i}}}& 0& 0& 0\\ 2\pi {\displaystyle \sum _{j=1, j\ne i}^N}{T}_{ij}& 0& 0& 0& 0& 0\\ {\displaystyle \frac{K_{{BESS}_i}}{T_{{BESS}_i}}}& 0& 0& 0& -{\displaystyle \frac{1}{T_{{BESS}_i}}}& 0\\ {\beta }_i& 0& 0& 1& 0& 0\end{array}\right]\\ \hfill {\mathcal{A}}_{ij}& =& \left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ -2\pi T_{ij}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\\ \hfill {\mathcal{B}}_i& =& {\left[\begin{array}{cccccc}0& 0& {\displaystyle \frac{1}{T_{g_i}}}& 0& 0& 0\end{array}\right]}^T\\ \hfill {\mathcal{C}}_i& =& \left[\begin{array}{cccccc}0& 0& 0& 0& 0& 1\\ {\beta }_i& 0& 0& 1& 0& 0\end{array}\right]\\ \hfill {\mathcal{D}}_i& =& \left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\end{array}\right]\\ \hfill {\mathcal{F}}_i& =& {\left[-\begin{array}{cccccc}{\displaystyle \frac{1}{T_{m_i}}}& 0& 0& 0& 0& 0\end{array}\right]}^T\end{array}$$

Considering that the data of measured outputs in each area are sending and receiving in the packet form, they are collected by digital devices, such as phasor measurement units (PMUs) and remote telemetry units (RTUs). Therefore, the discrete-time system model is discussed in this paper. The discrete-time representation of (10) can be expressed as

$$\left\{\begin{array}{l}{x}_i\left[k+1\right]=A_i{x}_i\left[k\right]+B_i{\mathsf{\Delta }P}_{c_i}\left[k\right]+{\displaystyle \sum _{j=1,j\ne i}^N}{A}_{ij}{x}_j\left[k\right]+F_i\left(\mathsf{\Delta }{P}_{{PV}_i}\left[k\right]+\mathsf{\Delta }{P}_{L_i}\left[k\right]\right)\\ y_i\left[k\right]=C_i{x}_i\left[k\right] \\ z_i\left[k\right]=D_i{x}_i\left[k\right]\end{array}\right.$$

(11)

where $A_i=e^{{\mathcal{A}}_i{T}_s}$, $B_i={\int }_0^{T_s}{e}^{{\mathcal{A}}_i\tau }{\mathcal{B}}_{i}d\tau$, $A_{ij}=e^{{\mathcal{A}}_{ij}{T}_s}$, $C_i={\mathcal{C}}_i$, $F_i={\int }_0^{T_s}{e}^{{\mathcal{A}}_i\tau }{\mathcal{F}}_{i}d\tau$, $D_i={\mathcal{D}}_i$, and $T_s$ is the sampling period.

3. Distributed LFC for Multi-Area Power Systems with Robust Performance

3.1. Proposed Controller

Considering that the state function combines the impacts of local dynamics and neighboring dynamics, the input of the controller also includes these two impacts. The controller of the ith area can be expressed as

$$\begin{array}{lll}{u}_i\left[k\right]& =& K_i{y}_i\left[k\right]+{\displaystyle \sum _{j=1, j\ne i}^N}{K}_{ij}{y}_j\left[k\right]\\ & =& {\mathsf{\Delta }P}_{c_i}\left[k\right]\end{array}$$

(12)

in which $K_i$ and $K_{ij}$ are controller gains for local area i and neighboring area j.

Thus, the derived closed-loop dynamics of area i is

$$\begin{array}{lll}{x}_i\left[k+1\right]& =& \left(A_i+B_i{K}_i{C}_i\right)x_i\left[k\right]\\ & +& \left({\displaystyle \sum _{j=1, j\ne i}^N}{A}_{ij}+B_i{K}_{ij}{C}_j\right)x_j\left[k\right]\\ & +& F_i(∆P_{{PV}_i}\left[k\right]+∆P_{L_i}[k\left]\right)\end{array}$$

(13)

3.2. $H_{\infty }$ Robust Performance Analysis

To design a robust LFC, the following definition needs to be first introduced.

Definition 1.

The multi-area power system (13) is asymptotically stable, for a given $H_{\infty }$ robust performance index $\gamma$, such that the following condition holds:

$${{‖∆f_i‖}_2}^2<{\gamma }^2{{‖∆P_{{total}_i}‖}_2}^2,=\forall ∆P_{{total}_i}\ne 0$$

(14)

where

$$∆P_{{total}_i}=∆P_{{PV}_i}+∆P_{L_i}$$

(15)

$${‖∆f_i‖}_2=\sqrt{{\displaystyle \sum _{k=0}^{\infty }}∆f_i^T\left[k\right]∆f_i\left[k\right]}$$

(16)

$${‖∆P_{{total}_i}‖}_2=\sqrt{{\displaystyle \sum _{k=0}^{\infty }}{∆P_{{total}_i}}^T\left[k\right]∆P_{{total}_i}\left[k\right]}$$

(17)

Then, Theorem 1 is introduced to show that the system in (13) is asymptotically stable for the frequency deviations $∆f_i$ under the $H_{\infty }$ robust performance index $\gamma$, if (14) is satisfied.

Theorem 1.

The multi-area power system (13) is asymptotically stable for the frequency deviations $∆f_i$ under the $H_{\infty }$ robust performance index $\gamma$, if there exist positive-definite matrices $M_i$ , i = 1, 2, …, N, such that

$$\left[\begin{array}{cc}{\tilde{A}}^{T}M\tilde{A}-M+D^{T}D& {\tilde{A}}^{T}MF\\ *& F^{T}MF-{\gamma }^{2}I\end{array}\right]<0$$

(18)

where * represents the symmetry term of a complex matrix,

$$\begin{array}{lll}\hfill x\left[k\right]& =& {\left[\begin{array}{ccc}{x}_1^T\left[k\right]& x_2^T\left[k\right]& \begin{array}{cc}\cdots & x_N^T\left[k\right]\end{array}\end{array}\right]}^T,\\ \hfill ∆P_L\left[k\right]& =& {\left[\begin{array}{ccc}∆P_{L_1}^T\left[k\right]& ∆P_{L_2}^T\left[k\right]& \begin{array}{cc}\cdots & ∆P_{L_N}^T\left[k\right]\end{array}\end{array}\right]}^T,\\ \hfill ∆f\left[k\right]& =& {\left[\begin{array}{ccc}∆f_1^T\left[k\right]& ∆f_2^T\left[k\right]& \begin{array}{cc}\cdots & ∆f_N^T\left[k\right]\end{array}\end{array}\right]}^T,\\ \hfill \tilde{A}& =& A+BKC,\\ \hfill B& =& diag\left\{\begin{array}{ccc}{B}_1& B_2& \begin{array}{cc}\cdots & B_N\end{array}\end{array}\right\},\\ \hfill D& =& diag\left\{\begin{array}{ccc}{D}_1& D_2& \begin{array}{cc}\cdots & D_N\end{array}\end{array}\right\},\\ \hfill F& =& diag\left\{\begin{array}{ccc}{F}_1& F_2& \begin{array}{cc}\cdots & F_N\end{array}\end{array}\right\},\\ \hfill M& =& diag\left\{\begin{array}{ccc}{M}_1& M_2& \begin{array}{cc}\cdots & M_N\end{array}\end{array}\right\},\\ \hfill A& =& \left[\begin{array}{cccc}{A}_1& A_{12}& \cdots & A_{1N}\\ A_{21}& A_2& \cdots & A_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ A_{N1}& A_{N2}& \cdots & A_N\end{array}\right],\\ \hfill K& =& \left[\begin{array}{cccc}{K}_1& K_{12}& \cdots & K_{1N}\\ K_{21}& K_2& \cdots & K_{2N}\\ ⋮& ⋮& \ddots & ⋮\\ K_{N1}& K_{N2}& \cdots & K_N\end{array}\right].\end{array}$$

(19)

Proof of Theorem 1.

According to the multi-area power system (13), define a Lyapunov function as follows:

$$V\left[k\right]=\sum _{i=1}^N{x}_i^T\left[k\right]M_i{x}_i\left[k\right]$$

(20)

where $M_i$ is the Lyapunov matrix.

In addition, according to (19), (13) can be written as

$$x\left[k+1\right]=\tilde{A}x\left[k\right]+F∆P_{total}\left[k\right]$$

(21)

Then, the forward difference of $V\left[k\right]$ is calculated as

$$\begin{array}{lll}∆V\left[k\right]& =& V\left[k+1\right]-V\left[k\right]\\ & =& {\displaystyle \sum _{i=1}^N}{x}_i^T\left[k+1\right]M_i{x}_i\left[k+1\right]-{\displaystyle \sum _{i=1}^N}{x}_i^T\left[k\right]M_i{x}_i\left[k\right]\end{array}$$

(22)

Combining (21) and (22), we can have

$$∆V\left[k\right]={\phi }^T\left[k\right]\left[\begin{array}{cc}{\tilde{A}}^{T}M\tilde{A}-M& {\tilde{A}}^{T}MF\\ *& F^{T}MF\end{array}\right]\phi \left[k\right]$$

(23)

where $\phi \left[k\right]={\left[\begin{array}{cc}{x}^T\left[k\right]& ∆{P_{total}}^T\left[k\right]\end{array}\right]}^T$.

When $∆P_{total}\left[k\right]=0$, we have

$$∆V\left[k\right]=x^T\left[k\right]{[\tilde{A}}^{T}M\tilde{A}-M\left]x\right[k]$$

(24)

Based on Theorem 1, it is obvious that $∆V\left[k\right]<0, \forall x\left[k\right]\ne 0$. Thus, the system (13) is asymptotically stable for $∆P_{total}\left[k\right]=0$.

When $∆P_{total}\left[k\right]\ne 0$, we design the performance criteria as

$$\psi \left[k\right]=V\left[k+1\right]-V\left[k\right]+∆f^T\left[k\right]∆f\left[k\right]-{\gamma }^2∆{P_{total}}^T\left[k\right]∆P_{total}\left[k\right]$$

(25)

Adding up both sides with $k=0, 1, \dots , \infty$, we have

$$\begin{array}{ll}\overline{\psi }={\displaystyle \sum _{k=0}^{\infty }}\psi \left[k\right]=& V\left[\infty \right]-V\left[0\right]+{\displaystyle \sum _{k=0}^{\infty }}∆f^T\left[k\right]∆f\left[k\right]\\ & -{\gamma }^2{\displaystyle \sum _{k=0}^{\infty }}∆{P_{total}}^T\left[k\right]∆P_{total}\left[k\right]\end{array}$$

(26)

According to (18), $\overline{\psi }<0$ can be derived. It is obvious that $V\left[\infty \right]>0$ and $V\left[0\right]=0$; therefore, ${{‖∆f_i‖}_2}^2<{\gamma }^2{{‖∆P_{{total}_i}‖}_2}^2$ always holds for any $∆P_{{total}_i}\ne 0$. Thus, Theorem 1 is proved. □

3.3. Design of the Proposed LFC

Though (18) indicates the conditions of the controller gains $K_i$ and $K_{ij}$, (18) is convex since $K_i$ and $K_{ij}$ are coupled with Lyapunov matrices $M_i$, which means the LMI tool in MATLAB cannot be employed directly. Thus, we proposed Theorem 2 to decouple the coupling in Theorem 1.

Theorem 2.

The multi-area power system (13) is asymptotically stable for the frequency deviations $∆f_i$ with the $H_{\infty }$ robust performance index $\gamma$, if there exist positive-definite matrices $M_i$, i = 1, 2, …, N, such that

$$\left[\begin{array}{ccc}-Z& \tilde{A}& F\\ *& D^{T}D-M& 0\\ *& *& -{\gamma }^{2}I\end{array}\right]<0$$

(27)

$$MZ=I$$

(28)

where $Z_i={M_i}^{-1}$.

Proof of Theorem 2.

Rewrite (18) as follows:

$$-\left[\begin{array}{c}{\tilde{A}}^T\\ F^T\end{array}\right]\left(-M\right)\left[\begin{array}{cc}\tilde{A}& F\end{array}\right]+\left[\begin{array}{cc}{D}^{T}D-M& 0\\ 0& -{\gamma }^{2}I\end{array}\right]<0$$

(29)

then employ the Schur complement, and we have

$$\left[\begin{array}{ccc}-M^{-1}& \tilde{A}& F\\ *& D^{T}D-M& 0\\ *& *& -{\gamma }^{2}I\end{array}\right]<0$$

(30)

Since $M$ is coupled with $M^{-1}$ in (28), letting $Z=M^{-1}$, and we have

$$\left[\begin{array}{ccc}-Z& \tilde{A}& F\\ *& D^{T}D-M& 0\\ *& *& -{\gamma }^{2}I\end{array}\right]<0$$

(31)

At this point, we can use LMI to obtain the controller gains $K_i$ and $K_{ij}$. □

4. Case Study

4.1. System Parameters

In the simulation, a four-area power system with PV and BESSs is adopted, and the interconnected topology is shown in Figure 2. Parameters of the thermal plant of four areas are listed in

Table 1. Parameter metrics.

Area 1	Area 2	Area 3	Area 4
$D_1=5$ pu/Hz	$D_2=1$ pu/Hz	$D_3=3$ pu/Hz	$D_4=4$ pu/Hz
$T_{m_1}=20$ pu·s	$T_{m_2}=14$ pu·s	$T_{m_3}=11$ pu·s	$T_{m_4}=9$ pu·s
$T_{{ch}_1}=1.2$ s	$T_{{ch}_2}=1.0$ s	$T_{{ch}_3}=0.7$ s	$T_{{ch}_4}=0.5$ s
$T_{g_1}=1.2$ s	$T_{g_2}=0.6$ s	$T_{g_3}=1.4$ s	$T_{g_4}=0.8$ s
$R_1=0.016$ Hz/pu	$R_2=0.03$ Hz/pu	$R_3=0.05$ Hz/pu	$R_4=0.04$ Hz/pu
$T_{12}=0.1$	$T_{21}=0.1$	$T_{31}=0.1$	$T_{41}=0.1$
$T_{13}=0.1$	$T_{23}=0.1$	$T_{32}=0.1$	$T_{42}=0$
$T_{14}=0.1$	$T_{24}=0$	$T_{34}=0.1$	$T_{43}=0.1$

[9].

Remark 1.

The transmission speed and transmission delay of real-time transmission are affected by the scalability of the power network. In practical power systems, it is usually divided into several areas, rather than dozens of areas. Thus, the chosen scalability of a four-area power system is reasonable in this paper.

4.2. Validation of Proposed DLFC for PV Variation

For PV power generation systems in four areas, we assume that the parameters of PV are the same, and they are chosen as follows [2]: $K_{{PV}_i}=0.2$, $T_{{PV}_i}=1.8$ s. The sampling time $T_s$ is set as 0.1 s, the $H_{\infty }$ robust performance index $\gamma$ is 1, and the initial condition $\mathsf{\Delta }f$ is 0.

The variation in solar radiation intensity $\mathsf{\Phi }$ will affect the PV output power and also cause fluctuation in the load frequency. In this case, the load disturbance $∆P_{L_i}$ is set to be zero to better demonstrate the ability of the proposed DLFC to mitigate the vibration of frequency deviations $∆f_i$. Assume that the locations of the four areas are the same and the solar radiation intensity received in the four areas is equal. The dynamics of frequency deviations with the variation in solar radiation intensity are shown in Figure 3, Figure 4 and Figure 5. $\mathsf{\Phi }$ in Figure 3 changes from 1 p.u. to 0.75 p.u. at t = 10 s, for Figure 4 it changes to 0.5 p.u., and for Figure 5 it changes to 0 p.u. We can notice that for different levels of variation in $\mathsf{\Phi }$, the proposed DLFC can restore load frequencies to the desired value when $\mathsf{\Phi }$ becomes stable, within small fluctuating limits.

4.3. Comparison Between Proposed Robust DLFC and Traditional PID

We assume that the parameters of the PV system and BESS are the same in the four areas, and the parameters are set as follows [2]: $K_{{PV}_i}=0.2$, $T_{{PV}_i}=1.8$ s, $K_{{BESS}_i}=-1/300$, and $T_{{BESS}_i}=0.1$ s. For the setting of robust DLFC in MATLAB R2024a, the $H_{\infty }$ robust performance index $\gamma$ is 1, and the load disturbance is $∆P_{L_i}=0.1\times \left(\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}\left(1\right)-0.5\right)$, in which rand(1) derives uniformly distributed random numbers with a range of [0, 1]. As for the sampling time $T_s$, a small sampling would allow the control system to accurately capture the dynamic behavior of the process, but after considering the computational efficiency, the sampling time $T_s=1 \mathrm{s}$ is chosen as a reasonable time to reach a balance between the need for precise control and the available computational resources.

Assume that solar radiation intensity received in the four areas is equal. Figure 6 illustrates the system dynamics of the PV system. The solar radiation is assumed to have several surges and drops in the simulation time. The solar radiation intensity $\mathsf{\Phi }$ is 1 p.u. during 0 s < t < 8 s, and then drops to 0.5 p.u. during 8 s < t < 20 s. It is found that $P_{PV}$ starts from 0 p.u. and increases to 0.2 p.u., and then decreases to 0.1 p.u. along with $\mathsf{\Phi }$. When 20 s < t < 40 s, $\mathsf{\Phi }$ stays at 0.5 p.u., and then rises to 1 p.u. during 40 s < t < 50 s. Same as $\mathsf{\Phi }$, $P_{PV}$ stays at 0.1 p.u. for 20 s and then increases to 0.2 p.u. After staying at 1 p.u. for 50 s < t < 70 s, $\mathsf{\Phi }$ uses 5 s to drop to 0.75 p.u., and stays there for 25 s. $P_{PV}$ also drops from 0.2 p.u. to 0.15 p.u. For 100 s < t < 110 s, $\mathsf{\Phi }$ keeps rising to 1 p.u. and finally stays at 1 p.u. until the simulation ends. $P_{PV}$ also increases to 0.2 p.u. The change in $P_{PV}$ is proportional to the variation in $\mathsf{\Phi }$.

Figure 7 shows the dynamics of frequency deviations $\mathsf{\Delta }{f}_i$ in each area. With the initial conditions set as $\mathsf{\Delta }{f}_i=0.02 \mathrm{H}\mathrm{z}$, the blue curve represents the proposed robust LFC, the red curve represents the conventional PID controller with controller gains $K_1=\left[-0.2585, 0.6189\right]$, $K_2=\left[-0.2568, 0.5947\right]$, $K_3=\left[-0.2594, 0.3945\right]$, $K_4=\left[-0.3859, 0.4851\right]$, and $K_{ij}=\left[0, 0\right]$ for $i,j=1,\dots , 4, i\ne j$, and the green curve represents no control in the power system.

Table 2. Distributed load frequency controller gains.

Area	Local Controller Gains	Neighboring Controller Gains
Area 1	K1 = [−0.1785, 0.5389]	K12 = [0.0194, −0.0085] K13 = [0.0295, −0.0171] K14 = [0.0266, −0.0036]
Area 2	K2 = [−0.2168, 0.5547]	K21 = [0.0089, 1.3977 × 10−4] K23 = [0.0229, −0.0072] K24 = [6.1100 × 10−4, −0.0027]
Area 3	K3 = [−0.1794, 0.3145]	K31 = [0.0070, 0.0013] K32 = [0.0138, −0.0023] K34 = [0.0121, 5.5528 × 10−4]
Area 4	K4 = [−0.3659, 0.2851]	K41 = [0.0104, 0.0016] K42 = [2.2252 × 10−4, −0.0026] K43 = [0.0220, 4.4162 × 10−4]

shows controller gains for the proposed robust DLFC. We can notice that with the proposed robust DLFC, the fluctuation in frequency deviation is lower than the conventional PID controller. Thus, with the proposed robust DLFC, the system is more stable and faster compared to the system with the conventional PID controller.

5. Sensitivity Analysis

5.1. The Impact of the Provided Capacity of PV Power

The proposed robust DLFC is sensitive to the penetration rate of the PV output power. The larger the provided capacity of PV power, the larger the fluctuation in frequency in the power system, and the more time needed to adjust. Increasing the capacity of PV power from 0.2 p.u. to 0.4 p.u. and 1 p.u., Figure 8, Figure 9 and Figure 10 show the corresponding dynamics of frequency deviations in each area. This case study was conducted without the impact of load disturbance, and the solar radiation intensity $\mathsf{\Phi }$ was the same as $\mathsf{\Phi }$ in the above case study. We can notice that as the capacity of PV power increases, the fluctuation in frequency deviation becomes larger, and it changes less smoothly. Though the capacity of PV power influences the performance of the controller, the proposed robust DLFC shows a better performance than the conventional PID controller, with a lower fluctuation in frequency deviation.

5.2. The Impact of Load Disturbance

In addition, the proposed robust DLFC is also sensitive to load disturbance. The large margin of the load disturbance and the high fluctuating frequency of the load disturbance will cause a larger fluctuation in frequency, and also a longer setting time to adjust the frequency dynamics. To explore the impact of different load disturbances on frequency change, in the following case study, the changing trend is kept and only the margin of load disturbance is modified, and the impact of PV output power is removed. Making the margin of load disturbance $∆P_{L_i}$ as 0.1, 0.15, and 0.2, Figure 11, Figure 12 and Figure 13 show the corresponding simulation results. We can notice that though the margin of load disturbance $∆P_{L_i}$ increases, the proposed robust DLFC always show a smaller fluctuation in frequency deviation compared to the conventional PID controller.

Figure 14 compares the dynamics of frequency deviations $\mathsf{\Delta }{f}_i$ in each area for different margins of load disturbance, as shown above. We can notice that with the increasing of the margin, the fluctuation in frequency deviations $\mathsf{\Delta }{f}_i$ is larger and frequency deviations $\mathsf{\Delta }{f}_i$ cross zero less frequently, which means that the controller needs more time to mitigate the fluctuation.

Furthermore, per Definition 1, with the $H_{\mathrm{\infty }}$ robust performance index $\gamma$, the power system is asymptotically stable when (14) is satisfied. By defining the $H_{\mathrm{\infty }}$ robust performance index $\gamma$, the impact of PV power perturbation and load disturbances on load frequency can be limited. According to (14), to reach a lower frequency deviation, the $H_{\mathrm{\infty }}$ robust performance index $\gamma$ may be defined as small to enhance robustness. However, if the chosen $H_{\mathrm{\infty }}$ robust performance index $\gamma$ is too small, the Lyapunov matrices $M_i$ may not exist, when selecting the controller gains by resorting to Theorem 2. Thus, the selection of the $H_{\mathrm{\infty }}$ robust performance index $\gamma$ needs to balance the solvability of the Lyapunov matrices $M_i$.

To quantitatively demonstrate the priority of the proposed robust DLFC over the traditional PID controller, we compare the mean square error (MSE) of the frequency deviation dynamic for margin $∆P_{L_i}=0.2$ under these two controllers. The formula of calculating MSE for $\mathsf{\Delta }{f}_i$ is specified as

$${MSE}_i={\displaystyle \frac{1}{S}}\sum _{k=1}^S{\mathsf{\Delta }{f}_i}^2\left(k\right)$$

(32)

where S indicates the quantity of the involved frequency deviation signal. Since the simulation time is 130 s and the sampling time is 1 s, S can be calculated as 130. We also calculate the percentage reduction for each area. The results are summarized in

Table 3. Comparison of MSE and percentage reduction in frequency deviation signal.

Area	MSE (Traditional Controller)	MSE (Proposed Controller)	Percentage Reduction
Area 1	2.9318 × 10−4	2.2358 × 10−5	92.37%
Area 2	6.9215 × 10−5	3.1932 × 10−5	53.87%
Area 3	1.3607 × 10−4	5.3123 × 10−5	60.96%
Area 4	6.9137 × 10−5	3.2527 × 10−5	52.95%

.

6. Conclusions

A robust DLFC scheme has been proposed for multi-area power systems with the integration of PV and BESSs. By introducing a performance index, two theorems have been formulated using LMIs to calculate proper controller gains with guaranteeing $H_{\infty }$ robust performance in the sense of asymptotic stability. Simulation results verify the robustness of the proposed controller against frequency fluctuation caused by the variation in PV, and demonstrate the effectiveness of the proposed method, compared with the traditional PID controller, by showing a lower fluctuation in frequency deviation than the traditional PID controller, and showing 92.37%, 53.87%, 60.96%, and 52.95% reductions under frequency deviation dynamics for a margin of $∆P_{L_i}=0.2$. In future work, we will consider incorporating some practical situations, such as transmission delay and cyber attacks, in multi-area power systems with PV and BESSs. Under these situations, it is more challenging for LFC to achieve robust performance, due to the more complicated modeling and stability analysis. To achieve this improvement, two contributions of this paper can be seen as the foundation.