Bi-layer optimal secondary frequency control approach for energy storage clusters considering wind power uncertainty

Increasing penetration of wind power with intermittency and variability threatens the stability of the power system frequency. The fast response capability of the energy storage system (ESS) makes it an effective measure to improve frequency regulation performance. However, designing an optimal control for numerous energy storage units (ESUs) with different power and energy characteristics is challenging. To solve the dilemma that the distributed control methods cannot achieve optimality over time horizons while the centralized optimization methods would cause high computational burdens, a bi-layer optimal control approach is proposed in this study. In the upper layer, a multi-area secondary frequency control (SFC) problem is established to determine control policies for ESS clusters under continuous wind power fluctuations and is solved by the Itô theory-based stochastic optimization (ITB-SO) method with high computational efficiency in a rolling-horizon manner. In the lower layer, the power output of ESUs is dispatched and coordinated using the distributed algorithm (DA) method via communication networks, considering different characteristics and the current state-of-charge (SoC) levels. Simulation studies carried out in a dual-region test system evidence the improvement of system frequency stability by compensating for rapid wind power fluctuations immediately via ESSs. The results show that the proposed ITB-SO method without scenario generation is suitable for real-time SFC for the ESS. The effectiveness of the DA method and its robustness in encountering loss of communication links are also verified.

1 Introduction

In the past few decades, the penetration of renewable energy sources (RES) has been increasing rapidly. At the end of 2020, the RES capacity, including wind power and photovoltaic power generation, was 1300GW, accounting for 9.4% of the total installed power generation capacity in the world (Khan et al., 2021). The increasing RES penetration level in the interconnected power grids results in a decrease in the system inertia (Meegahapola and Flynn, 2010). For example, the wind turbine generators are decoupled from the system frequency by converters, so they cannot provide frequency regulation support for a long time (Tsili and Papathanassiou, 2009). Severer frequency fluctuations would happen in the areas with low system inertia than others in the interconnected power systems (Miller et al., 2011). Meanwhile, the intermittency and variability nature of RES power generation would lead to continuous power imbalances and system frequency deviations (Bidram and Davoudi, 2012). Therefore, the increasing integration of RES will inevitably bring challenges to the stability of system frequency.

Exploiting energy storage systems (ESSs) for maintaining frequency stability of the systems with high penetration levels of the RES has recently attracted attention both in industry and academia. In recent years, the application of the ESS is gradually becoming widespread. In 2021, the rated power of installed ESS was 173.7GW, and 1,363 ESS projects were operational globally (U.S.DOE, 2021). The ESS has higher ramping and faster response capability than conventional generators, so better dynamic tracking performance is available for frequency regulation control. Since the traditional automatic generation control (AGC) signal is unsuitable for ESSs with limited energy capacity, filtering mechanisms are developed to extract fast-moving and high-frequency components of the AGC signal for the ESS to participate in SFC (Cheng et al., 2014). Meanwhile, the management of state-of-charge (SoC) for ESS during frequency regulation is important as well. It needs to prevent the second frequency drop issue due to sudden withdrawals of some ESSs. Karrari et al. (2020) propose an SoC recovery mechanism for flywheel ESS so that the residual energy capacity is always in the optimal state. Garcia-Torres et al. (2021) apply the stochastic optimization (SO) method considering the management SoC level and address the prediction uncertainties of RES.

In order to further exploit the complementary advantages of different types of ESSs in terms of rated power and energy capacity, hybrid ESSs are developed and utilized in frequency regulation. The main ESS types, including batteries (Li-ion, Pb-Acid), flywheels, and supercapacitors, possess different power-to-energy ratio characteristics. Jan et al. (2020) used a fuzzy PI controller for ESSs to cater to the frequency variations adaptively. Shim et al. (2018) proposed droop control with the SoC feedback to improve their frequency response and regulation services to the grid. However, the control strategies in these studies are separate and independent, which means different ESSs cannot be coordinated, and hence unnecessary mutual regulation and waste of resources would be caused. Esmaili et al. (2013) decomposed the power signal into fast/slow components to allocate responsibility for hybrid ESSs. Oshnoei et al. (2020) adopted the model predictive control method to determine an ESS control scheme with various characteristics, which is robust to uncertain disturbances. Nevertheless, these centralized ways are sensitive to noise and communication failure (Zhao and Ding, 2018b). On the contrary, the distributed realizations relying on local and neighboring information can reduce computation time and take advantage of local intelligent agents (Anderson et al., 2021).

Distributed methods such as the consensus algorithm are widely used in economic dispatch problems (Li et al., 2019), control for photovoltaic arrays (Zhang et al., 2021) and batteries (Chen C. et al., 2019), etc. Alsharif et al. (2020) applied the distributed algorithm (DA) method coordinating ESS power output to enhance the frequency nadir and improve the dynamic performance. Lee et al. (2016) coordinated the group battery ESS by a distributed control algorithm for voltage and frequency deviation regulation. Cherukuri and Cortes (2018) developed a distributed method for ESSs considering the change of load, but uniform storage efficiency is assumed. However, these studies did not take into account random wind power generation or the influence of ESS control strategies on the frequency of interconnected grids. Zhao and Ding (2018a) present a DA method to maximize the total welfare of battery ESSs during frequency regulation considering changes in wind power. But, it relies on accurate predictions of wind power, which is difficult to accomplish, considering their inherent variability. As a result, the stochastic characteristics of RES should be taken into account when designing accurate control strategies for the ESS.

However, a disadvantage of the distribution methods is that the optimality of power dispatching between ESSs can only be achieved at the current time step, instead of considering the optimal solution over a future finite time horizon. The power output scheme is possibly not adapted to the energy-constrained ESSs when the signal is constantly biased, resulting in storage saturation. It indicates that the ESS can participate in the SFC only if the energy content of the signal is kept low (Megel et al., 2018). So, it is necessary to develop reasonable control strategies to maintain optimal energy levels of ESSs. But, on the other hand, if the power output of all ESSs is decided in a centralized optimization problem over a period, a large scale of variables would lead to computational burdens and intractability, which is unsuitable for real-time SFCs.

In order to fill the abovementioned research gap, a bi-layer SFC framework for the ESS is presented in this article. Different energy storage units (ESUs) in the same area are combined as an ESS cluster to reduce the computational complexities of upper-layer control, and the optimal control policies are updated in a rolling-horizon manner to ensure frequency regulation performance. Then, the ESUs coordinate their power output via communication networks in real-time. Therefore, different from the previous studies, the proposed approach can achieve a balance between the optimization method considering wind power uncertainty and the distribution method with high efficiency. The main contributions of this article are summarized as follows:

1) An upper-layer optimization problem of multi-area frequency regulation considering stochastic wind power fluctuation is established to determine the optimal control policy of ESS clusters. The Itô theory-based stochastic optimization (ITB-SO) method is proposed to efficiently obtain the solutions, in which the wind power uncertainty is modeled by stochastic differential equations (SDE) without scenario generation. The computational burden is significantly reduced so that it is suitable for real-time frequency regulation.

2) To make full use of the local computing capability, a DA method with high efficiency and robustness is applied to dispatching power output among ESUs. The power command can be adaptively decomposed considering the different characteristics of power-type and energy-type units. The evaluation of the SoC level is improved by adopting the logistic function to better reflect the safe and dangerous SoC zones. Sudden withdrawal of ESUs due to energy saturation and violation of power limits can be avoided to guarantee power tracking performance.

The remainder of the article is organized as follows: Section 2 describes the multi-area frequency regulation model considering ESSs and wind power uncertainties. Section 3 proposes an optimal SFC approach for ESS clusters based on the ITB-SO method. Section 4 proposes a power dispatching approach for ESUs in the clusters based on the DA method. Case studies are presented in Section 5. Section 6 gives a conclusion.

2 System model

2.1 Dynamic model of frequency regulation considering the ESS

The system frequency is affected by load and generation sources, including conventional generators and ESSs. Regardless of the power loss, the dynamic frequency response model of multi-areas can be expressed as follows (Mauricio et al., 2009; Li et al., 2016):

$\frac{d}{dt}{\Delta f}_i=-\frac{\mathrm{1}}{\mathrm{2}{H}_i}\left(D_i{\Delta f}_i-{\Delta p}_i^{\mathrm{M}}+p_i^{\mathrm{E}}+{\Delta p}_i^{\mathrm{L}}+{\displaystyle \sum }_{j:i\to j}{\Delta p}_{ij}^{\mathrm{t}\mathrm{i}\mathrm{e}}\right),\forall i,j\in \mathcal{A}(1)$

$\frac{d}{dt}{\Delta p}_i^{\mathrm{M}}=-\frac{\mathrm{1}}{T_i^{\mathrm{t}}}\left({\Delta p}_i^{\mathrm{M}}-{\Delta p}_i^{\mathrm{A}\mathrm{G}\mathrm{C}}+\frac{\mathrm{1}}{R_i}{\Delta f}_i\right),\forall i\in \mathcal{A}(2)$

$\frac{d}{dt}{\Delta p}_{ij}^{\mathrm{t}\mathrm{i}\mathrm{e}}=2\pi \sum T_{ij}\left({\Delta f}_i-{\Delta f}_j\right),\forall i,j\in \mathcal{A}(3)$

$\frac{d}{dt}{\Delta p}_i^{\mathrm{A}\mathrm{G}\mathrm{C}}=-K_i^{\mathrm{I}}\left(B_i{\Delta f}_i+{\displaystyle \sum }_{j:i\to j}{\Delta p}_{ij}^{\mathrm{t}\mathrm{i}\mathrm{e}}\right),\forall i,j\in \mathcal{A}(4)$

where $\mathcal{A}$ is the set of area systems; ${\Delta f}_i$ is the frequency deviation of the ith area from the rated value; ${\Delta p}_i^{\mathrm{M}}$ is the active power deviation of conventional generators from the nominal value in the ith area; $p_i^{\mathrm{E}}$ is the active power of the ith ESS cluster because numerous units can be equivalently combined to reduce the number of variables; ${\Delta p}_i^{\mathrm{L}}$ is the change of load; ${\Delta p}_{ij}^{\mathrm{t}\mathrm{i}\mathrm{e}}$ represents the tie line power from the ith to the jth area; ${\Delta p}_i^{\mathrm{A}\mathrm{G}\mathrm{C}}$ is the AGC reference; $H_i$ and $D_i$ are the equivalent system inertia and damping, respectively; $T_i^{\mathrm{t}}$ is turbine time constant; $R_i$ is the droop coefficient; $T_{ij}$ is the synchronizing coefficient; $B_i$ is the frequency bias factor; and $K_i^{\mathrm{I}}$ is the integral coefficient. Eq 1 gives the frequency dynamic characteristic; Eq. 2 shows the change of regulation output from the generators, including primary frequency regulation and AGC; Eq. 3 shows the power change on the tie line between different areas; tie-line bias control in Eq. 4 automatically adjusts the AGC reference by driving the area control errors (ACE) to zero.

The wind power output can be considered a negative system load as shown in Eq. (5), which is more difficult to accurately forecast than load. Hence, the continuous wind power fluctuations would pose a threat to the safety and stability of the system frequency.

${\Delta p}_i^{\mathrm{L}}=-{\Delta p}_i^{\mathrm{W}},\forall i\in \mathcal{A}(5)$

where ${\Delta p}_i^{\mathrm{W}}$ is the deviation of wind power from its initial value in the ith area.

2.2 Energy storage model

Different ESUs are combined into an ESS cluster, so the overall power output and energy level of the ith cluster can be represented by the sum of that of the member ESUs as follows:

$p_i^{\mathrm{E}}=p_i^{\mathrm{E},\mathrm{d}}-p_i^{\mathrm{E},\mathrm{c}}={\displaystyle \sum }_k\left(p_k^{\mathrm{d}}-p_k^{\mathrm{c}}\right),\forall k\in {\mathcal{S}}_i(6)$

$s_i^{\mathrm{E}}={\displaystyle \sum }_k{s}_k,\forall k\in {\mathcal{S}}_i(7)$

$\frac{d}{dt}{s}_k=p_k^{\mathrm{c}}{\eta }_k^{\mathrm{c}}-p_k^{\mathrm{d}}/{\eta }_k^{\mathrm{d}},\left({\eta }_k^{\mathrm{c}},{\eta }_k^{\mathrm{d}}\right)\in (\mathit{0},\mathit{1}],\forall k\in {\mathcal{S}}_i(8)$

where $\mathcal{S}$ is the set of ESS clusters; $p_i^{\mathrm{E},\mathrm{d}}$ and $p_i^{\mathrm{E},\mathrm{c}}$ are the discharge and charge power of the ith ESS cluster, respectively; $p_k^{\mathrm{d}}$ and $p_k^{\mathrm{c}}$ are the discharge and charge power of the kth ESU, respectively; $s_i^{\mathrm{E}}$ and $s_k$ are the energy level of the ESS cluster and ESU, respectively; the initial condition is $s_{\mathit{0}}$; and ${\eta }_k^{\mathrm{c}}$ and ${\eta }_k^{\mathrm{d}}$ are the charge and discharge efficiencies, respectively. The power output and energy level satisfy the limit constraints as shown in (9 and 10), and the state-of-charge (SoC) can be obtained according to the energy level as shown in (11).

$0\le p_k^{\mathrm{d}}\le P_k^{\max },0\le p_k^{\mathrm{c}}\le P_k^{\max }(9)$

$\mathit{0}\le s_k\le S_k^{\max }(10)$

${SoC}_k=\frac{s_k}{S_k^{\max }}(11)$

where $P_k^{\max }$ represents the rated power; $S_k^{\max }$ is the rated energy capacity; and ${SoC}_k$ is the SoC level of the kth ESU.

2.3 Model of wind power prediction error

SDE has been widely used as an effective model of RES prediction error, such as wind power (Verdejo et al., 2016) and photovoltaic power (Lingohr and Müller, 2019). The probability distribution information and the temporal correlation characteristic could be simultaneously expressed. It is convenient to fit the non-Gaussian distribution via the Itô process by adopting different functions in SDE models so as to well-reflect the detailed characteristics of RES uncertainty. On the other hand, the SDE model can be embedded into SO problems without scenario generation in order to significantly release computational burdens.

Specifically, the wind power output value consists of the ultra-short-term prediction in a finite time period and the prediction error, which is modeled via the Itô process as follows:

${\Delta p}_i^{\mathrm{W}}={\Delta p}_i^{\mathrm{W},\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}+{\xi }_i,\forall i\in \mathcal{A}(12)$

$d{\xi }_i=\mu \left({\xi }_i\right)dt+\sigma \left({\xi }_i\right)d{W}_t,\forall i\in \mathcal{A}(13)$

$\mu \left({\xi }_i\right)=-\frac{\left({\xi }_i-b_{\mathit{1}}\right)}{\tau },\sigma \left({\xi }_i\right)=\frac{b_{\mathit{2}}}{\tau }\left({\xi }_i-b_{\mathit{3}}\right)\left(b_{\mathit{4}}-{\xi }_i\right)(14)$

where ${\Delta p}_i^{\mathrm{W},\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}$ is the predicted wind power deviating from the initial forecast; ${\xi }_i$ represents the prediction error; $\mu \left(\cdot \right)$ and $\sigma \left(\cdot \right)$ are the drift and diffusion function of SDE, respectively; $W_t$ is a Wiener process; $\tau$ represents the time constant; parameter $b_{\mathit{1}}$ is the mean deviation from the predicted value; $b_{\mathit{2}}$ is the fluctuation intensity; and $b_{\mathit{3}}$ and $b_{\mathit{4}}$ indicate the fluctuation interval. The parameters can be easily adjusted in different time periods according to the preference of operators in advance.

2.4 Bi-layer SFC framework for ESS clusters

The optimal SFC for ESS clusters is divided into two control layers. In the upper layer, the convex SO problems are established to decide the optimal power output minimizing the frequency deviation. However, modeling numerous ESUs will result in a large number of decision variables for the optimization problem, which is difficult to solve rapidly. Therefore, optimal control policies for the overall power output of ESS clusters are determined to reduce computational complexities. The ITB-SO method is applied instead of scenario-based methods in order to quickly obtain the optimal solutions for the SFC.

In the lower layer, the total power command is dispatched between the ESUs in the same ESS cluster with a communication network. The optimal dispatching scheme is obtained by the DA method considering the current SoC levels and different characteristics of power-type and energy-type units. In summary, the framework of the bi-layer optimal SFC is shown in Figure 1.

FIGURE 1

FIGURE 1. Framework of bi-layer secondary frequency control for energy storage clusters in multi-area systems.

3 Optimal SFC of ESSs based on ITB-SO

3.1 Objective function

The frequency regulation control problem to minimize the frequency deviation and operation cost of the ESS can be described as a SO problem over $\mathcal{T}:t\in \left[\mathit{0},T\right]$, and the objective function is defined as

$\min J={\mathbb{E}}_{{\mathit{\xi }}_{\mathit{0}}}\left\{{\int }_{\mathrm{0}}^T\left({\Delta \mathit{f}}_t^{\mathrm{T}}{\mathbf{\Lambda }}_f{\Delta \mathit{f}}_t+{\mathit{R}}_t^{\mathrm{T}}{\mathbf{\Lambda }}_R{\mathit{R}}_t\right)dt\right\}+{\mathbb{E}}_{{\mathit{\xi }}_{\mathit{0}}}\left\{{\mathrm{\Lambda }}_S{\displaystyle \sum }_i{\left({SoC}_{i,T}-{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right)}^2\right\}(15)$

where $\mathcal{T}$ is the set of time index; T represents the terminal time; ${\Delta \mathit{f}}_t$ is the vector form of the frequency deviation of each area shown in (16); ${\mathit{R}}_t$ is the vector of control variables of the ESS shown in (17), including the discharge and charge components of power output. Since the objective is to achieve a minimum of the square power output, simultaneous non-zero charging and discharging power can be avoided; ${SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}$ is the reference value for SoC, which is set as 0.5 in this study; ${\mathit{\xi }}_t$ is the vector of wind power prediction error ${\xi }_i$; ${\mathbb{E}}_{{\mathit{\xi }}_{\mathit{0}}}\left\{\cdot \right\}$ denotes the expectation operator under the initial conditions ${\mathit{\xi }}_{\mathit{0}}$; ${\mathbf{\Lambda }}_f$, ${\mathbf{\Lambda }}_R$, and ${\Lambda }_S$ are the weight parameters.

${\Delta \mathit{f}}_t={\left[{\Delta f}_{i,t}\right]}^{\mathrm{T}},\forall i\in \mathcal{A},\forall t\in \mathcal{T}(16)$

${\mathit{R}}_t={\left[p_{i,t}^{\mathrm{E},\mathrm{d}},p_{i,t}^{\mathrm{E},\mathrm{c}}\right]}^{\mathrm{T}},\forall i\in \mathcal{A},\forall t\in \mathcal{T}(17)$

Since different ESUs are combined into an ESS cluster and meet the power command, their overall SoC level can be estimated by

$\frac{d}{dt}{SoC}_i=\frac{1}{S_i^{\max }}\left(p_i^{\mathrm{E},\mathrm{c}}{\eta }_i^{\mathrm{c},\mathrm{e}\mathrm{q}}-p_i^{\mathrm{E},\mathrm{d}}/{\eta }_i^{\mathrm{d},\mathrm{e}\mathrm{q}}\right),\forall i\in \mathcal{A}(18)$

where ${\eta }_i^{\mathrm{d},\mathrm{e}\mathrm{q}}$ and ${\eta }_i^{\mathrm{c},\mathrm{e}\mathrm{q}}$ are the equivalent discharge and charge efficiency, respectively, which are set as the average over ESUs; and $S_i^{\max }$ is the total rated energy capacity of the ith ESS cluster.

It can be seen that the objective function (15) shows the trade-off between frequency regulation performance and the power and energy cost of the ESS. At the terminal time, the SoC level is forced to be close to the reference value by the objective function because it is necessary to maintain the feasibility of the next rolling horizon after the terminal time.

Moreover, the frequency deviation constraint and SoC constraint should be considered, as shown in (19) and (20), to enhance frequency regulation performance and guarantee the operation of the ESS. It is to be noted that constraint (19) could be violated sometimes when ESSs do not have enough power capacity to participate in frequency regulation, and hence it should be removed if the problem is infeasible.

$-{\Delta f}^{\max }\le {\Delta f}_{i,t}\le {\Delta f}^{\max },\forall i\in \mathcal{A},\forall t\in \mathcal{T}(19)$

$\underset{¯}{SoC}\le {SoC}_{i,t}\le \overline{SoC},\forall i\in \mathcal{A},\forall t\in \mathcal{T}(20)$

3.2 Itô theory-based stochastic optimization method

3.2.1 Series expansion of stochastic assessment function

To obtain the expected value of the objective function under all uncertainties, an ITB approximation method is adopted to evaluate it in a deterministic way (Chen et al., 2019c) so that the SO problem can be solved efficiently. For convenience, the stochastic system is expressed as (21) considering the stochastic resources ${\mathit{\xi }}_t$ described by SDE.

$\begin{array}{c}d{\mathit{X}}_t=\left(\mathit{A}{\mathit{X}}_t+{\mathit{B}\mathit{\xi }}_t\right)dt\\ d{\mathit{\xi }}_t=\mu \left({\mathit{\xi }}_t\right)dt+\sigma \left({\mathit{\xi }}_t\right)d{W}_t\end{array}(21)$

where ${\mathit{X}}_t$ represents the state variables, which is a $N_x\times \mathit{1}$ column vector; ${\mathit{\xi }}_t$ is a $N_{\xi }\times \mathit{1}$ column vector representing wind power prediction error in this study; A and B are the coefficient matrices; and the initial condition of the stochastic system is ${\mathit{X}}_{\mathit{0}}$ and ${\mathit{\xi }}_0$.

The stochastic assessment function (SAF) is defined and denoted as $u\left(t,{\mathit{X}}_{\mathit{0}}\right)$ when specific functions of ${\mathit{X}}_t$ are given and can be expanded into a series of deterministic terms as follows:

$u\left(t,{\mathit{X}}_{\mathit{0}}\right)={\mathbb{E}}_{{\mathit{\xi }}_0}\left\{{\int }_0^{t}g\left({\mathit{X}}_s\right)ds+h\left({\mathit{X}}_t\right)\right\}=\underset{n}{{\displaystyle \sum }^{\infty }}{\stackrel{\sim }{u}}_n\left(t,{\mathit{X}}_{\mathrm{0}}\right)(22)$

where $g\left(\cdot \right)$ and $h\left(\cdot \right)$ are functions of ${\mathit{X}}_t$ with continuous second-order derivatives and ${\stackrel{\sim }{u}}_n$ is the nth order deterministic components of the SAF. Therefore, SAF can be approximately estimated by ignoring the high-order terms.

Specifically, the objective function (15) evaluated by SAF at time T can be decomposed into low-order terms as

$J=u\left(T,{\mathit{X}}_{\mathrm{0}}\right)\approx {\stackrel{\sim }{J}}_{\mathrm{0}}+{\stackrel{\sim }{J}}_{\mathrm{1}}(23)$

Given that $g={\Delta \mathit{f}}_t^{\mathrm{T}}{\mathbf{\Lambda }}_f{\Delta \mathit{f}}_t+{\mathit{R}}_t^{\mathrm{T}}{\mathbf{\Lambda }}_R{\mathit{R}}_t$ and $h={\mathrm{\Lambda }}_S{\displaystyle \sum }_i{\left({SoC}_{i,t}-{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right)}^2$, the 0-th and 1-st order terms of the objective function are defined as

${\stackrel{\sim }{J}}_{\mathit{0}}={\int }_0^{T}g\left({\stackrel{\sim }{\mathit{X}}}_t\right)dt+h\left({\stackrel{\sim }{\mathit{X}}}_T\right)(24)$

${\stackrel{\sim }{J}}_{\mathit{1}}={\int }_0^T\frac{\mathit{1}}{\mathit{2}}{\sigma \left({\stackrel{\sim }{\mathit{\xi }}}_t\right)}^{\mathrm{T}}{\mathit{\theta }}_t\sigma \left({\stackrel{\sim }{\mathit{\xi }}}_t\right)dt(25)$

where ${\stackrel{\sim }{\mathit{X}}}_t$ is the auxiliary state variables by ignoring $\sigma \left({\mathit{\xi }}_t\right)d{W}_t$ in (21). Hence, they can be recognized as the expected state variables in the deterministic system without the influence of uncertain resources; ${\stackrel{\sim }{\mathit{\xi }}}_t$ is independent of the state variables, so $\sigma \left({\stackrel{\sim }{\mathit{\xi }}}_t\right)$ can be obtained once the initial condition ${\mathit{\xi }}_0$ is given. ${\mathit{\theta }}_t$ represents the second-order derivative of SAF components and is defined as

${\mathit{\theta }}_t={\int }_{\mathrm{0}}^t{\widehat{\mathit{X}}}_s^{\mathrm{T}}\left[{\nabla }^{\mathrm{2}}g\right]{\widehat{\mathit{X}}}_{s}ds+{\widehat{\mathit{X}}}_t^{\mathrm{T}}\left[{\nabla }^{\mathrm{2}}h\right]{\widehat{\mathit{X}}}_t(26)$

where ${\nabla }^{\mathit{2}}$ is the operator of the Hessian matrices and ${\widehat{\mathit{X}}}_t$ is another group of auxiliary state variables with ${\widehat{\mathit{X}}}_t=\partial {\stackrel{\sim }{\mathit{X}}}_t/\partial {\mathit{\xi }}_0,{\widehat{\mathit{X}}}_{\mathit{0}}={\mathit{0}}_{N_{\xi }\times N_x}$. The detailed proof of the abovementioned series expansion method based on Itô theory can be found in (Chen et al., 2019c).

3.2.2 Control policy

The control output of the ESS should be appropriately adjusted in response to different wind power values. However, it is hard to accurately predict wind power so the design of the control strategy is challenging considering unrevealed uncertainties. Disturbance feedback control is needed to address the problem (Skaf and Boyd, 2010). An effective measure is parameterizing the control policies as an affine policy (Chen et al., 2019b; Qiu et al., 2020), which is expressed in (27). The control policy only needs to be updated once in each rolling period, and once ${\mathit{\xi }}_t$ is measured, the control output can be determined in real-time.

${\mathit{R}}_t=r\left({\mathit{\xi }}_t\right)={\mathit{r}}_{\mathit{0},t}+\mathit{K}{\mathit{\xi }}_t,\forall t\in \mathcal{T}(27)$

where ${\mathit{r}}_{\mathit{0},t}$ is the base control output and K is the gain matrix which leads to correction of power output according to the disturbance ${\mathit{\xi }}_t$. Thus, the auxiliary control variables are

$\begin{array}{c}{\stackrel{\sim }{\mathit{R}}}_t={\mathit{r}}_{\mathit{0},t}+\mathit{K}{\stackrel{\sim }{\mathit{\xi }}}_t,\forall t\in \mathcal{T}\\ {\widehat{\mathit{R}}}_t=\mathit{K}{\widehat{\mathit{\xi }}}_t,\forall t\in \mathcal{T}\end{array}(28)$

where ${\stackrel{\sim }{\mathit{R}}}_t$ and ${\widehat{\mathit{R}}}_t$ are included in ${\stackrel{\sim }{\mathit{X}}}_t$ and ${\widehat{\mathit{X}}}_t$, respectively; there is ${\widehat{\mathit{\xi }}}_t=\partial {\stackrel{\sim }{\mathit{\xi }}}_t/\partial {\mathit{\xi }}_0,{\widehat{\mathit{\xi }}}_{\mathit{0}}={\mathit{I}}_{N_{\xi }\times N_{\xi }}$. Based on the control variables, the auxiliary variables of SoC can also be obtained as

$\begin{array}{c}\frac{d}{dt}{\stackrel{\sim }{\mathit{S}\mathit{o}\mathit{C}}}_t=\frac{1}{{\mathit{S}}^{\max }}\mathit{\eta }{\stackrel{\sim }{\mathit{R}}}_t,\forall t\in \mathcal{T}\\ \frac{d}{dt}{\widehat{\mathit{S}\mathit{o}\mathit{C}}}_t=\frac{1}{{\mathit{S}}^{\max }}\mathit{\eta }\mathit{K}{\widehat{\mathit{\xi }}}_t,\forall t\in \mathcal{T}\end{array}(29)$

where ${\mathit{S}\mathit{o}\mathit{C}}_t$ and $\mathit{\eta }$ are the vector forms of ${SoC}_{i,t}$ and $\left[{\eta }_i^{\mathrm{c},\mathrm{e}\mathrm{q}},\mathit{1}/{\eta }_i^{\mathrm{d},\mathrm{e}\mathrm{q}}\right]$, respectively. It can be seen that the auxiliary variable of SoC is related to the historical information of wind power uncertainties, so the assessment of SoC would be gradually inaccurate over time. Therefore, the rolling-horizon manner should be used to guarantee good control performance.

3.2.3 Chance constraint

In the established SO problem, the constraints are classified into equality constraints, including (1)–(6), (18), (28), and (29); and inequality constraints, including (9), (19), and (20). The linear inequality constraints are considered in a probabilistic way as chance constraints, expressed as

$\mathrm{Pr}\left\{{\mathit{\varphi }}_r^{\mathrm{T}}{\mathit{X}}_t\le {\overline{\varphi }}_r\right\}\ge \gamma ,\forall r\in \mathcal{C}(30)$

where $\mathcal{C}$ is the set of all inequality constraints; $\mathrm{Pr}\left\{\cdot \right\}$ represents the probability operator; ${\mathit{\varphi }}_r$ and ${\overline{\varphi }}_r$ are the coefficients vector and the upper bound of the rth constraint, respectively; and $\gamma$ is the chance tolerance. Chance constraint (30) can be rewritten as (31) for inner approximation with coefficient ${\kappa }_{\gamma }$ under specific probability level (Calafiore and Ghaoui, 2006).

${\kappa }_{\gamma }\sqrt{\mathrm{v}\mathrm{a}\mathrm{r}\left\{{\mathit{\varphi }}_r^{\mathrm{T}}{\mathit{X}}_t\right\}}+{\mathit{\varphi }}_r^{\mathrm{T}}{\stackrel{\sim }{\mathit{X}}}_t\le {\overline{\varphi }}_r(31)$

According to the abovementioned series expansion theory, the variance term can be evaluated by SAF as (32), and a detailed explanation of (31) and (32) can be found in Supplementary Appendix SA.

$\mathrm{v}\mathrm{a}\mathrm{r}\left\{{\mathit{\varphi }}_r^{\mathrm{T}}{\mathit{X}}_t\right\}={\int }_0^t{\sigma \left({\stackrel{\sim }{\mathit{\xi }}}_s\right)}^{\mathrm{T}}{\widehat{\mathit{X}}}_s^{\mathrm{T}}{\mathit{\varphi }}_r{\mathit{\varphi }}_r^{\mathrm{T}}{\widehat{\mathit{X}}}_s\sigma \left({\stackrel{\sim }{\mathit{\xi }}}_s\right)ds(32)$

3.2.4 Convex optimization problem reformulation

The convexity of the SO problem can be guaranteed when the affine feedback policy is adopted, and the optimization problem is formed as a quadratically constrained quadratic programming (QCQP) problem. The trapezoidal rule (Sanchez-Gasca et al., 1995) can be used in the numerical calculation of the integral. To express the chance constraints as a second-order cone constraint form, the coefficient vector ${\mathit{\pi }}_r$ is introduced to reformulate the variance term in (32), considering the accumulation characteristic of integral operator as

$\begin{array}{c}\sqrt{{\mathit{\pi }}_r^{\mathrm{T}}{\mathit{Y}}_t{\mathit{Y}}_t^{\mathrm{T}}{\mathit{\pi }}_r}\le -\frac{\mathit{1}}{{\kappa }_{\gamma }}\left({\mathit{\varphi }}_r^{⊺}{\stackrel{\sim }{\mathit{X}}}_t-{\overline{\varphi }}_r\right),\forall r\in \mathcal{C}\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}{\mathit{Y}}_t={\left[{\widehat{\mathit{X}}}_{\mathit{0}}^{\mathrm{T}},{\widehat{\mathit{X}}}_{\mathit{1}}^{\mathrm{T}},\cdots ,{\widehat{\mathit{X}}}_t^{\mathrm{T}}\right]}^{\mathrm{T}},\forall t\in \mathcal{T}\end{array}(33)$

Moreover, the ‘ = ’ in the objective function can be replaced by ‘≥’ to relax and can be rewritten into second-order cone constraints similar to the method in (33). In summary, the overall optimization problem is expressed as

$\begin{array}{c}\underset{{\mathit{r}}_{\mathit{0},t},\mathit{K},{\stackrel{\sim }{\mathit{X}}}_t,{\widehat{\mathit{X}}}_t}{\min }J\ge {\phi }^{\mathit{2}}\\ \begin{array}{cc}\mathrm{s}.\mathrm{t}.& \begin{array}{c}{\parallel {\mathit{\omega }}_{\mathit{0}}^{\mathrm{T}}\stackrel{\sim }{\mathit{X}}\parallel }_{\mathrm{2}}\le {\phi }_{\mathit{0}},{\parallel {\mathit{\omega }}_{\mathit{1}}^{\mathrm{T}}\widehat{\mathit{X}}\parallel }_{\mathrm{2}}\le {\phi }_{\mathit{1}},{||\begin{array}{c}{\phi }_{\mathit{0}}\\ {\phi }_{\mathit{1}}\end{array}||}_{\mathrm{2}}\le \phi \\ {\parallel {\mathit{\pi }}_r^{\mathrm{T}}{\mathit{Y}}_t\parallel }_{\mathrm{2}}\le -\frac{\mathit{1}}{{\kappa }_{\gamma }}\left({\mathit{\varphi }}_r^{⊺}{\stackrel{\sim }{\mathit{X}}}_t-{\overline{\varphi }}_r\right),\forall r\in \mathcal{C}\\ \mathrm{e}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\left(1\right)-\left(6\right),\left(18\right),\left(28\right),\left(29\right)\end{array}\end{array}\end{array}(34)$

where ${\phi }_{\mathit{0}}$ and ${\phi }_{\mathit{1}}$ are the square root value of ${\stackrel{\sim }{J}}_{\mathit{0}}$ and ${\stackrel{\sim }{J}}_{\mathit{1}}$, respectively; ${\mathit{\omega }}_{\mathit{0}}$ and ${\mathit{\omega }}_{\mathit{1}}$ are the introduced coefficient vectors; and ${\parallel \cdot \parallel }_{\mathit{2}}$ represents 2-norm operator.

4 Power dispatching approach for ESUs based on DA

Once the total power command is given by the upper-layer optimal SFC, different ESUs in the same ESS cluster are required to coordinate with each other to determine their power output and completely meet the command. However, the characteristics of different ESUs are probably inconsistent, including rated power, rated energy capacity, and the current SoC level. Thus, the possible energy saturation and violation of power limits would lead to sudden withdrawals of ESUs, which might influence the frequency regulation effect. Therefore, a DA method is proposed in this section to coordinate the power output of ESUs with communication networks and guarantee the frequency regulation performance.

4.1 Optimization control problem of ESUs

Similar to the cost of conventional generators in economic dispatch (Wen et al., 2018), the virtual cost function of ESUs can be defined as a quadratic function with respect to the power output of the ESU and the current SoC level (Megel et al., 2018) at each time step t as

$C_{k,t}\left(p_{k,t}\right)=\frac{\mathrm{1}}{\mathrm{2}}{\alpha }_k{p}_{k,t}^{\mathit{2}}+{\beta }_k\left({SoC}_{k,t}\right)p_{k,t}(35)$

where $p_{k,t}$ denotes $p_{k,t}^{\mathrm{d}}-p_{k,t}^{\mathrm{c}}$; ${\alpha }_k$ is the quadratic term coefficient; and ${\beta }_k\left({SoC}_{k,t}\right)$ is the linear term coefficient, used as an evaluation function of the current SoC level. The detailed derivation of (35) can be found in Supplementary Appendix SB.

Thus, the optimization control problem of ESUs is to achieve the minimizing objective function (36) under equality constraint (37) and other inequality constraints.

$\min \underset{k=\mathrm{1}}{{\displaystyle \sum }^{N_i}}{C}_{k,t}\left(p_{k,t}\right)(36)$

$\underset{k=\mathrm{1}}{{\displaystyle \sum }^{N_i}}{p}_{k,t}=p_{i,t}^{\mathrm{E}},\forall i\in \mathcal{A},\forall k\in {\mathcal{S}}_i(37)$

where $N_i$ is the number of ESUs in the ith ESS cluster and $p_{i,t}^{\mathrm{E}}$ is the total power command determined by the upper control layer for all ESUs in the same cluster.

In addition, the inequality constraints are twofold: power and energy limits, and power direction constraints, i.e., $\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(p_{k,t}\right)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(p_{i,t}^{\mathrm{E}}\right)$. It is to avoid the opposite direction of the ESU power output, which might cause some to charge others and result in a waste of resources.

This problem is a quadratic programming problem and can be solved in a centralized way (Yin et al., 2021). However, the local computing capability in a parallel way is not fully used, so the efficiency of solving optimization problems would be reduced. Therefore, the DA needs to be developed.

According to the general process of utilizing distributed methods to solve optimization problems (Wang et al., 2019), the Lagrange function $L\left(\mathit{p},\lambda \right)$ is constructed in (38) based on the objective function (36) and equality constraint (37). $\mathit{p}$ represents the vector form of $p_k$. The inequality constraints are ignored here, which will be considered later in DA. For convenience, the time step label t is omitted hereafter.

$L\left(\mathit{p},\lambda \right)=\underset{k=\mathrm{1}}{{\displaystyle \sum }^{N_i}}{C}_k\left(p_k\right)+\lambda \left(p_i^{\mathrm{E}}-\underset{k=\mathrm{1}}{{\displaystyle \sum }^{N_i}}{p}_k\right)(38)$

where $\lambda$ is the Lagrange multiplier. If $\partial L\left(\mathit{p},\lambda \right)/\partial p_k=\mathit{0},\forall k$, for any ESU, there is

$\frac{\partial C_k\left(p_k\right)}{\partial p_k}-\lambda =\mathit{0},\forall k(39)$

In other words, when all $\partial C_k\left(p_k\right)/\partial p_k$ are equal to $\lambda$, the original objective function (36) can achieve its extreme value. Therefore, the value of the consensus variable $x_k$ should approach the value of $\lambda$ as shown in (40). In addition, $x_k$ is the input for the dispatching scheme of $p_k$, where a higher $x_k$ leads to a higher $p_k$ and vice versa.

$x_k=\frac{\partial C_k\left(p_k\right)}{\partial p_k}={\alpha }_k{p}_k+{\beta }_k\left({SoC}_k\right),\forall k(40)$

Specifically, the evaluation function ${\beta }_k\left({SoC}_k\right)$ with the basic form of logistic function (Postnikov, 2020) is proposed in this article to better-reflect the safe and dangerous zones of SoC as shown in (41).

${\beta }_k\left({SoC}_k\right)=w_k\cdot \mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left({SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}-{SoC}_k\right)\cdot F\left({SoC}_k\right)(41)$

where $w_k$ is the weight factor with $w_k\in \left[\mathit{0,1}\right]$ and $F\left({SoC}_k\right)$ is constructed above and below the reference value and defined as

$F\left({SoC}_k\right)=\frac{e^{m\left(\left|{SoC}_k-{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right|-\frac{\mathrm{1}}{\mathrm{2}}\left|{SoC}_k^{\lim }-{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right|\right)}}{\mathrm{1}+e^{m\left(\left|{SoC}_k-{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right|-\frac{\mathrm{1}}{\mathrm{2}}\left|{SoC}_k^{\lim }-{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}\right|\right)}}(42)$

where ${SoC}_k^{\lim }$ is the SoC limit: when ${SoC}_k<{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}$, there is ${SoC}_k^{\lim }={SoC}_k^{\min }$; when ${SoC}_k>{SoC}^{\mathrm{r}\mathrm{e}\mathrm{f}}$, there is ${SoC}_k^{\lim }={SoC}_k^{\max }$; and m is an adaptive factor used to adjust the steepness of $F\left({SoC}_k\right)$ in (42), in order to adjust the scope of safe and dangerous SoC levels. After several comparative analyses, m is selected as 35 as shown in Figure 2.

FIGURE 2

FIGURE 2. Different curves of F (SoC) when the value of the adaptive factor m increases from 25 to 65.

4.2 Discrete consensus algorithm

The main idea of DA is that the agents use a communication network to share information to complete the iterative calculation. The communication network topology can be modeled by a directed or undirected graph denoted by $G=\left(V,E\right)$ (Pourbabak et al., 2018), where $V$ is the set of nodes and $E$ is the set of edges formed by adjacent nodes. A Laplace matrix $\mathit{L}=\left[L_{kn}\right]$ is defined as shown in (43), where the diagonal element is the degree of nodes in graph $G$, and the non-diagonal element is the (0–1) element $a_{kn}$, which is determined by whether the nodes k and n are adjacent.

$\{\begin{array}{c}{L}_{kk}={\displaystyle \sum }_{k\ne n}{a}_{kn}\\ L_{kn}=-a_{kn}\end{array}(43)$

The agents obtain information from their neighbors through the communication network and update their own information. All agents finally converge to the same value after several iterations. The first-order discrete consensus algorithm can be updated according to (44).

$x_k\left[l+\mathrm{1}\right]=\underset{n=\mathrm{1}}{{\displaystyle \sum }^N}{d}_{kn}\left[l\right]x_n\left[l\right](44)$

where l is the discrete sequence representing the lth iteration; $x_k$ is the consensus variable of the kth agent; and $d_{kn}\left[l\right]$ is the element of ${\mathit{D}}_{n\times n}={[d}_{kn}]$ and is defined as

$d_{kn}\left[l\right]=\frac{\left|L_{kn}\right|}{\underset{n=\mathrm{1}}{{\displaystyle \sum }^N}\left|L_{kn}\right|}(45)$

4.3 Control procedure of DA

1) Step 1: Initializing power dispatch and consensus variables.

Set $l=\mathit{0}$. Current SoC level information is measured to compute the evaluation function by (41). The total power command is allocated equally to each ESU as the initial condition, and the consensus variable is initialized by (40).

2) Step 2: Iterative calculation.

The consensus variables are iterated by DA. In order to ensure that the final dispatching scheme can match the total power command, a power correction term is introduced for the leader agent. The update rule is shown as

$\{\begin{array}{c}{x}_k\left[l+\mathrm{1}\right]=\underset{n=\mathrm{1}}{{\displaystyle \sum }^{N_i}}{d}_{kn}\left[l\right]x_n\left[l\right]+\delta p_{\mathrm{d}\mathrm{e}\mathrm{v}}\left[l\right],k\mathrm{i}\mathrm{s}\mathrm{l}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{e}\mathrm{r}\\ x_k\left[l+\mathrm{1}\right]=\underset{n=\mathrm{1}}{{\displaystyle \sum }^{N_i}}{d}_{kn}\left[l\right]x_n\left[l\right],k\mathrm{i}\mathrm{s}\mathrm{f}\mathrm{o}\mathrm{l}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}\end{array}(46)$

where $\delta$ is the convergence coefficient that determines the convergence speed of the leader and $p_{\mathrm{d}\mathrm{e}\mathrm{v}}$ is the deviation from command at the current iteration:

$p_{\mathrm{d}\mathrm{e}\mathrm{v}}\left[l\right]=p_i^{\mathrm{E}}-\underset{n=\mathrm{1}}{{\displaystyle \sum }^{N_i}}{p}_k\left[l\right](47)$

$p_k$ is updated according to the consensus variable $x_k$ by the following rules:

${\overline{p}}_k\left[l\right]=\frac{x_k\left[l\right]-{\beta }_k\left({SoC}_k\right)}{{\alpha }_k}(48)$

$p_k\left[l\right]={\left\{{\overline{p}}_k\left[l\right]\right\}}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}(49)$

where (49) makes a restriction on the theoretical power output by (48), which is possibly infeasible due to inequality constraints, i.e., power and energy limits and power direction constraints.

Then, the remaining power deviation would be fed back to the consensus variable through the correction term and shared by other neighbors.

3) Step 3: Judging the convergence condition.

First, each ESU agent should reach an agreement with the adjacent agents on their consensus variables. However, complete consensus cannot be achieved for all agents actually, so the approximate condition is implemented as

$\left|x_k\left[l+\mathrm{1}\right]-x_n\left[l+\mathrm{1}\right]\right|<{\epsilon }_{\mathrm{1}},k\leftrightarrow n(50)$

where $k\leftrightarrow n$ means adjacent and ${\epsilon }_{\mathit{1}}$ is a small constant.

Second, the power deviation is also taken as a factor in the convergence condition as shown in (51). In order to avoid over-iteration, ${\epsilon }_{\mathit{2}}$ is set as a small constant close to 0, and it determines the accuracy of tracking the command.

$\left|p_{\mathrm{d}\mathrm{e}\mathrm{v}}\left[l+\mathrm{1}\right]\right|<{\epsilon }_{\mathrm{2}}(51)$

In summary, the flowchart of the control procedure of the DA method is shown in Figure 3.

FIGURE 3

FIGURE 3. Flowchart of the DA method.

5 Case study

5.1 Simulation model setting

The simulations are performed using Matlab software and run on a machine equipped with Intel i7 3.2 GHz CPU and 16 GB RAM. The established optimization problems are solved by the Yalmip toolbox called Gurobi optimizer. The test simulations based on the actual wind power scenario are conducted using the Simulink platform.

The structure of the simulation system for frequency regulation is shown in Figure 4 containing two areas with rated capacities of both 600 MW and a rated system frequency of both 50 Hz. The parameters of the simulation system are listed in Table 1 (Kundur et al., 1994; Pathak et al., 2018). The parameters and initial SoC levels of the ESSs are listed in Table 2. The predicted and actual wind power deviations in the two areas are shown in Figure 5. The total simulation time is set as 7200s; the time step is set as 1s; the parameters of the SDE models can be found in Supplementary Appendix SC; the upper-layer SO problem is solved every 60 s, in which the parameter ${\Delta f}^{\max }$ is set as 0.2Hz, ${\kappa }_{\gamma }$ as 1.7, and $\underset{¯}{SoC}$ and $\overline{SoC}$ as 0.2 and 0.8, respectively; and the lower-layer DA is implemented at each time step, in which the matrices ${\mathit{D}}_i$ are listed in Supplementary Appendix SD.

FIGURE 4

FIGURE 4. Structure of two-area systems considering wind power and energy storage systems.

TABLE 1

TABLE 1. Parameters of the system frequency regulation model in p.u.

TABLE 2

TABLE 2. Parameters of energy storage units.

FIGURE 5

FIGURE 5. Predicted and actual wind power deviation in (A) area 1 and (B) area 2.

5.2 Frequency regulation results

In this section, two contrast cases are introduced for comparison. Contrast method 1: the contribution of the ESS is not considered during the frequency regulation, and Contrast method 2, namely, PI: the power command for the ESS is computed via PI controller, where the input signal of PI is high-frequency components of the ACE signal via a high-pass filter. In addition, our proposed method is denoted as ITB-SO for convenience.

The results of frequency deviations under these different control methods are compared in Figure 6. It is clear that when the ESS participates in frequency regulation, the peak and nadir value of system frequency is always reduced and increased. Most of the time, the frequency deviations of the two areas can be limited within the ±0.2 Hz bounds, both under the PI and the ITB-SO methods. However, the PI method cannot respond immediately to wind power fluctuation, so it is difficult to make full use of the rapid adjustment capability of the ESS. The root mean squared (RMS) value of frequency deviations is given for comparison as shown in Figure 6. It can be seen that the RMS value under the ITB-SO method is the lowest, and that of the two different areas can be restricted to a similar level.

FIGURE 6

FIGURE 6. Comparison of frequency deviation under different control methods of (A) area 1 and (B) area 2.

To verify that the proposed ITB-SO method can adaptively adjust the power output of the ESS considering different rated power and energy capacity, Figure 7 gives the simulation results under the PI and the ITB-SO methods. It can be seen that the power command for ESS cluster 2 cannot be accurately met by the actual output during several periods because the rated power of ESS cluster 2 is small. On the contrary, the scheme determined by ITB-SO is based on the estimation of the future possible power and SoC in a finite horizon; therefore, the violation of power and energy restriction can be avoided, and the frequency regulation performance can be guaranteed.

FIGURE 7

FIGURE 7. Comparison results under different control methods of (A) ESS cluster 1 and (B) ESS cluster 2.

In addition, the SoC level under ITB-SO is closer to the reference value than that of the PI method. For example, the ESS cluster 2 operates near the lower SoC bound, i.e., 0.2 after 5500 s under PI control, while under ITB-SO, ESS cluster 1, with a higher energy level, discharges more power to release the burden of ESS cluster 2. That is because the penalty of SoC deviation is also considered in the objective function in order to force the SoC to the reference value.

5.3 Accuracy and results of stochastic assessment

The effectiveness of the rolling-horizon manner is validated in this section by comparing it with Monte Carlo sample (MCS) simulations, where the MCS simulations are conducted under 100 wind power generation scenarios. Results of the two successive rollings (100th and 101st) are shown in Figure 8. The violet bands represent 85–100% (5% interval) of all 100 MCS results, and the actual operation results are clearly in the scope of the MCS. The four subgraphs in the left column show the wind power deviation, frequency deviation, power, and SoC of the ESS cluster in area 1, and those in area 2 are in the right column.

FIGURE 8

FIGURE 8. Actual results and Monte Carlo samples of (A,B) wind power deviation, (C,D) frequency deviation (E,F) power output of the ESS cluster, and (G,H) SoC level of the ESS cluster.

In Figures 8C and D, the expected value is obtained without taking into account wind power uncertainty, where the impact of stochastic fluctuation on the objective and chance constraints in the optimization problem cannot be considered. Therefore, the application of SO is necessary. In Figures 8E and F, it can be seen that the control output of the ESS can be further corrected according to the feedback policy and respond to the real-time wind power prediction error rapidly. In Figures 8G and H, the curve marked ‘real’ is the mean value of all ESUs in the ESS cluster, and the estimated results are obtained by equivalent discharge and charge efficiency as in (18). It indicates that accuracy can be ensured, and the effectiveness of ITB-SO can be verified.

To verify the effectiveness of the approximately decomposed reformulation of SAF in the first-order form, the objective value assessed by SAF of each rolling iteration is compared with the results via MCS simulations. The gross objective is divided into multi-objective for comparison as shown in Figure 9, including the objective of frequency deviations, control variables, and SoC deviations. It evidenced the accuracy of SAF without scenario generation; thus, the computational efficiency can be enhanced.

FIGURE 9

FIGURE 9. Objective value estimated by SAF compared with Monte Carlo samples: (A) objective of frequency deviations, (B) objective of control variables, and (C) objective of SoC deviations. Keywords: energy storage, frequency regulation, wind power uncertainty, stochastic optimization, Itô theory, distributed algorithm.

5.4 Results of DA

Taking the ESUs in the ESS cluster 1 as an example, Figure 10 shows their output power and SoC level under the DA method, compared with the contrast method, which is to reach the average command signals. It is easy to be implemented without iteration, and the power command can be met under ideal situations.

FIGURE 10

FIGURE 10. Power output and SoC level of each ESU compared under (A,B) DA control and (C,D) equally allocation method.

However, the tracking error cannot be eliminated because of the SoC saturation or violation of rated power in some cases when the power command is equally allocated to each ESU. For example, in Figures 10C and D, some ESUs withdraw from frequency regulation since their SoC level reaches the upper bound (such as 800 s) and the lower bound (such as 2500 s) or the power signal is larger than the rated output (such as 1080 s), while others with large energy capacity are barely used. Therefore, the effect of frequency regulation could deteriorate due to sudden withdrawals. In contrast, the proposed DA method enables all ESUs to coordinate their output power and adjust their SoC at optimal levels. It is clear from Figure 10B that the SoC level is attracted to their reference value and the saturation can be avoided. Therefore, the DA method can make full use of the characteristics and advantages of power-type and energy-type units effectively.

Random communication failure between the ESU agents is implemented to verify the robustness of the DA dispatching method. The results are studied in different simulation systems containing 5, 10, 20, and 50 ESUs. The situation with complete connectivity is defined as all units being connected with each other. Then, the proportion of missing communication links varies from 10 to 80% (10% interval) by disconnecting the communication links randomly. The boxplot of iteration times and the mean value of algebraic connectivity are shown in Figure 11. The algebraic connectivity is the second smallest eigenvalue of the Laplace matrix $\mathit{L}$, which reflects the connectivity of the graph (Yoonsoo and Mesbahi, 2006).

FIGURE 11

FIGURE 11. Results of iteration times and mean value of connectivity for simulation systems with different numbers of ESUs Nomenclature.

When parts of the communication links disconnect between ESUs, they can still share information through other remaining links and successfully track the power command as well. It can be found in Figure 11 that with the increase in the number of missing links, the connectivity gradually decreases while iteration times increase. However, the trend is diverse with different numbers of units: iteration times with 20 and 50 units increase slower than that with 5 and 10 units, and the iteration characteristic is more stable. It indicates that random communication failure has a slight impact on the system with numerous ESUs and the robustness is higher.

6 Conclusion

In this article, a bi-layer optimal SFC approach is proposed for ESSs to participate in multi-area systems frequency regulation under continuous wind power fluctuations. The upper-layer control based on the ITB-SO method in a rolling-horizon manner is implemented to determine the power command for different ESS clusters considering wind power uncertainties. The power output of ESUs in the same cluster is coordinated using the DA method in the lower-layer control to track the total power command. The simulation applications reveal that 1) the stability of the system frequency can be enhanced by compensating for wind power volatility immediately via ESSs. Moreover, the trade-off between frequency regulation performance and the operation cost of the ESS can be achieved. 2) The ITB-SO method can solve optimization problems precisely without generating scenarios of wind power. The computational efficiency is high, so it is suitable for real-time SFCs. 3) Sudden withdrawals of ESSs due to energy saturation and violation of power limit can be avoided by dispatching power signals using the DA method. The robustness of the DA method is also validated while missing limited communication links between ESUs. In future work, the addition of PV systems would be considered.

Data availability statement

The original contributions presented in the study are included in the article; further inquiries can be directed to the corresponding author.

Author contributions

ZW, FW, and LS contributed to the conception and design of the study. ZW and YL are responsible for providing experimental design and data analysis. ZW and WH are responsible for providing code implementation. All authors contributed to the manuscript writing revised, read, and approved the submitted version.

Funding

This work was supported in part by the Funds for International Cooperation and Exchanges of the National Natural Science Foundation of China under Grant 52061635102, and in part by the Overseas Expertise Introduction Project for Discipline Innovation (“111Project”) under Grant B14022.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg.2022.1005281/full#supplementary-material

References

Alsharif, H., Jalili, M., and Hasan, K. N. 2020. "Analyzing the impact of distribution of battery energy storage system for participation in frequency regulation", in: 2020 Australasian Universities Power Engineering Conference (AUPEC), Hobart, Australia, 29 November 2020 - 02 December 2020, 1–6.

Google Scholar

Anderson, T., Muralidharan, M., Srivastava, P., Haghi, H. V., Cortés, J., Kleissl, J., et al. (2021). Frequency regulation with heterogeneous energy resources: A realization using distributed control. IEEE Trans. Smart Grid 12 (5), 4126–4136. doi:10.1109/TSG.2021.3071778

CrossRef Full Text | Google Scholar

Bidram, A., and Davoudi, A. (2012). Hierarchical structure of microgrids control system. IEEE Trans. Smart Grid 3 (4), 1963–1976. doi:10.1109/TSG.2012.2197425

CrossRef Full Text | Google Scholar

Calafiore, G. C., and Ghaoui, L. E. (2006). On distributionally robust chance-constrained linear programs. J. Optim. Theory Appl. 130 (1), 1–22. doi:10.1007/s10957-006-9084-x

CrossRef Full Text | Google Scholar

Chen, C., Bao, Y.-Q., Wu, X.-H., Wang, B., and Shen, C. (2019a). Battery energy storage system based on incremental cost consensus algorithm for the frequency control. IEEE Access 7, 147362–147372. doi:10.1109/access.2019.2946458

CrossRef Full Text | Google Scholar

Chen, X., Lin, J., Liu, F., and Song, Y. (2019b). Optimal control of AGC systems considering non-Gaussian wind power uncertainty. IEEE Trans. Power Syst. 34 (4), 2730–2743. doi:10.1109/tpwrs.2019.2893512

CrossRef Full Text | Google Scholar

Chen, X., Lin, J., Liu, F., and Song, Y. (2019c). Stochastic assessment of AGC systems under non-Gaussian uncertainty. IEEE Trans. Power Syst. 34 (1), 705–717. doi:10.1109/tpwrs.2018.2865502

CrossRef Full Text | Google Scholar

Cheng, Y., Tabrizi, M., Sahni, M., Povedano, A., and Nichols, D. (2014). Dynamic available AGC based approach for enhancing utility scale energy storage performance. IEEE Trans. Smart Grid 5 (2), 1070–1078. doi:10.1109/TSG.2013.2289380

CrossRef Full Text | Google Scholar

Cherukuri, A., and Cortes, J. (2018). Distributed coordination of DERs with storage for dynamic economic dispatch. IEEE Trans. Autom. Contr. 63 (3), 835–842. doi:10.1109/tac.2017.2731809

CrossRef Full Text | Google Scholar

Esmaili, A., Novakovic, B., Nasiri, A., and Abdel-Baqi, O. (2013). A hybrid system of Li-ion capacitors and flow battery for dynamic wind energy support. IEEE Trans. Ind. Appl. 49 (4), 1649–1657. doi:10.1109/tia.2013.2255112

CrossRef Full Text | Google Scholar

Garcia-Torres, F., Bordons, C., Tobajas, J., Real-Calvo, R., Santiago, I., and Grieu, S. (2021). Stochastic optimization of microgrids with hybrid energy storage systems for grid flexibility services considering energy forecast uncertainties. IEEE Trans. Power Syst. 36 (6), 5537–5547. doi:10.1109/tpwrs.2021.3071867

CrossRef Full Text | Google Scholar

Jan, M. U., Xin, A., Abdelbaky, M. A., Rehman, H. U., and Iqbal, S. (2020). Adaptive and fuzzy PI controllers design for frequency regulation of isolated microgrid integrated with electric vehicles. IEEE Access 8, 87621–87632. doi:10.1109/access.2020.2993178

CrossRef Full Text | Google Scholar

Karrari, S., Baghaee, H. R., Carne, G. D., Noe, M., and Geisbuesch, J. (2020). Adaptive inertia emulation control for high‐speed flywheel energy storage systems. IET Gener. Transm. Distrib. 14 (22), 5047–5059. doi:10.1049/iet-gtd.2020.0066

CrossRef Full Text | Google Scholar

Khan, F. H., Pal, T., Kundu, B., and Roy, R. "Wind energy: A practical power analysis approach", in: 2021 Innovations in Energy Management and Renewable Resources,Kolkata, 07 February 2021, 1–6.

CrossRef Full Text | Google Scholar

Kundur, P., Balu, N. J., and Lauby, M. G. (1994). Power system stability and control. New York: McGraw-Hill.

Google Scholar

Lee, S.-J., Kim, J.-H., Kim, C.-H., Kim, S.-K., Kim, E.-S., Kim, D.-U., et al. (2016). Coordinated control algorithm for distributed battery energy storage systems for mitigating voltage and frequency deviations. IEEE Trans. Smart Grid 7 (3), 1713–1722. doi:10.1109/tsg.2015.2429919

CrossRef Full Text | Google Scholar

Li, N., Zhao, C., and Chen, L. (2016). Connecting automatic generation control and economic dispatch from an optimization view. IEEE Trans. Control Netw. Syst. 3 (3), 254–264. doi:10.1109/tcns.2015.2459451

CrossRef Full Text | Google Scholar

Li, Q., Gao, D. W., Zhang, H., Wu, Z., and Wang, F. y. (2019). Consensus-based distributed economic dispatch control method in power systems. IEEE Trans. Smart Grid 10 (1), 941–954. doi:10.1109/TSG.2017.2756041

CrossRef Full Text | Google Scholar

Lingohr, D., and Müller, G. (2019). Stochastic modeling of intraday photovoltaic power generation. Energy Econ. 81, 175–186. doi:10.1016/j.eneco.2019.03.007

CrossRef Full Text | Google Scholar

Mauricio, J. M., Marano, A., Gomez-Exposito, A., and Martinez Ramos, J. L. (2009). Frequency regulation contribution through variable-speed wind energy conversion systems. IEEE Trans. Power Syst. 24 (1), 173–180. doi:10.1109/tpwrs.2008.2009398

CrossRef Full Text | Google Scholar

Meegahapola, L., and Flynn, D. 2010. "Impact on transient and frequency stability for a power system at very high wind penetration", in: IEEE PES General Meeting, Minneapolis, 29 July 2010, 1–8.

CrossRef Full Text | Google Scholar

Megel, O., Liu, T., Hill, D. J., and Andersson, G. (2018). Distributed secondary frequency control algorithm considering storage efficiency. IEEE Trans. Smart Grid 9 (6), 6214–6228. doi:10.1109/tsg.2017.2706979

CrossRef Full Text | Google Scholar

Miller, N. W., Clark, K., and Shao, M. "Frequency responsive wind plant controls: Impacts on grid performance", in: 2011 IEEE Power and Energy Society General Meeting, Detroit, 28 July 2011, 1–8.

Google Scholar

Oshnoei, A., Kheradmandi, M., and Muyeen, S. M. (2020). Robust control scheme for distributed battery energy storage systems in load frequency control. IEEE Trans. Power Syst. 35 (6), 4781–4791. doi:10.1109/tpwrs.2020.2997950

CrossRef Full Text | Google Scholar

Pathak, N., Bhatti, T. S., Verma, A., and Nasiruddin, I. (2018). AGC of two area power system based on different power output control strategies of thermal power generation. IEEE Trans. Power Syst. 33 (2), 2040–2052. doi:10.1109/tpwrs.2017.2734923

CrossRef Full Text | Google Scholar

Postnikov, E. B. (2020). Estimation of COVID-19 dynamics “on a back-of-envelope”: Does the simplest SIR model provide quantitative parameters and predictions? Chaos Solit. Fractals 135, 109841. doi:10.1016/j.chaos.2020.109841

PubMed Abstract | CrossRef Full Text | Google Scholar

Pourbabak, H., Luo, J., Chen, T., and Su, W. (2018). A novel consensus-based distributed algorithm for economic dispatch based on local estimation of power mismatch. IEEE Trans. Smart Grid 9 (6), 5930–5942. doi:10.1109/tsg.2017.2699084

CrossRef Full Text | Google Scholar

Qiu, Y., Lin, J., Liu, F., Song, Y., Chen, G., and Ding, L. (2020). Stochastic online generation control of cascaded run-of-the-river hydropower for mitigating solar power volatility. IEEE Trans. Power Syst. 35 (6), 4709–4722. doi:10.1109/tpwrs.2020.2991229

CrossRef Full Text | Google Scholar

Sanchez-Gasca, J. J., D'Aquila, R., Price, W. W., and Paserba, J. J. 1995. "Variable time step, implicit integration for extended-term power system dynamic simulation", in: Power Industry Computer Application Conference, Salt Lake City, 2 May 1995.

Google Scholar

Shim, J. W., Verbic, G., Zhang, N., and Hur, K. (2018). Harmonious integration of faster-acting energy storage systems into frequency control reserves in power grid with high renewable generation. IEEE Trans. Power Syst. 33 (6), 6193–6205. doi:10.1109/tpwrs.2018.2836157

CrossRef Full Text | Google Scholar

Skaf, J., and Boyd, S. P. (2010). Design of affine controllers via convex optimization. IEEE Trans. Autom. Contr. 55 (11), 2476–2487. doi:10.1109/TAC.2010.2046053

CrossRef Full Text | Google Scholar

Tsili, M., and Papathanassiou, S. (2009). A review of grid code technical requirements for wind farms. IET Renew. Power Gener. 3 (3), 308–332. doi:10.1049/iet-rpg.2008.0070

CrossRef Full Text | Google Scholar

U.S.DOE (2021). Global energy storage database projects. [Online]. Available: https://sandia.gov/ess-ssl/gesdb/public/projects.html (accessed 01Jun, 2022).

Google Scholar

Verdejo, H., Awerkin, A., Saavedra, E., Kliemann, W., and Vargas, L. (2016). Stochastic modeling to represent wind power generation and demand in electric power system based on real data. Appl. Energy 173, 283–295. doi:10.1016/j.apenergy.2016.04.004

CrossRef Full Text | Google Scholar

Wang, Y., Xu, Y., Tang, Y., Liao, K., Syed, M. H., Guillo-Sansano, E., et al. (2019). Aggregated energy storage for power system frequency control: A finite-time consensus approach. IEEE Trans. Smart Grid 10 (4), 3675–3686. doi:10.1109/TSG.2018.2833877

CrossRef Full Text | Google Scholar

Wen, G., Yu, X., Liu, Z.-W., and Yu, W. (2018). Adaptive consensus-based robust strategy for economic dispatch of Smart grids subject to communication uncertainties. IEEE Trans. Ind. Inf. 14 (6), 2484–2496. doi:10.1109/tii.2017.2772088

CrossRef Full Text | Google Scholar

Yin, J., Leon, J. I., Perez, M. A., Franquelo, L. G., Marquez, A., and Vazquez, S. (2021). Model predictive control of modular multilevel converters using quadratic programming. IEEE Trans. Power Electron. 36 (6), 7012–7025. doi:10.1109/tpel.2020.3034294

CrossRef Full Text | Google Scholar

Yoonsoo, K., and Mesbahi, M. (2006). On maximizing the second smallest eigenvalue of a state-dependent graph Laplacian. IEEE Trans. Autom. Contr. 51 (1), 116–120. doi:10.1109/tac.2005.861710

CrossRef Full Text | Google Scholar

Zhang, X., Yu, T., Yang, B., and Jiang, L. (2021). A random forest-assisted fast distributed auction-based algorithm for hierarchical coordinated power control in a large-scale PV power plant. IEEE Trans. Sustain. Energy 12 (4), 2471–2481. doi:10.1109/tste.2021.3101520

CrossRef Full Text | Google Scholar

Zhao, T., and Ding, Z. (2018a). Cooperative optimal control of battery energy storage system under wind uncertainties in a microgrid. IEEE Trans. Power Syst. 33 (2), 2292–2300. doi:10.1109/tpwrs.2017.2741672

CrossRef Full Text | Google Scholar

Zhao, T., and Ding, Z. (2018b). Distributed agent consensus-based optimal resource management for microgrids. IEEE Trans. Sustain. Energy 9 (1), 443–452. doi:10.1109/TSTE.2017.2740833

CrossRef Full Text | Google Scholar

Glossary

RES Renewable energy sources

SFC Secondary frequency control

ESS Energy storage system

ESU Energy storage unit

SoC State-of-charge

AGC Automatic generation control

ACE Area control errors

ITB-SO Itô theory-based stochastic optimization

DA Distributed algorithm

SDE Stochastic differential equation

SAF Stochastic assessment function

MCS Monte Carlo samples

Variables

${\Delta f}_i$ Frequency deviation of the ith area

${\Delta p}_i^{\mathrm{M}}$ Active power deviation of conventional generators in the ith area

$p_i^{\mathrm{E}}$ Active power of ith ESS cluster

${\Delta p}_i^{\mathrm{L}}$ Change of load in the ith area

${\Delta p}_{ij}^{\mathrm{t}\mathrm{i}\mathrm{e}}$ Tie line power from the ith to the jth area

${\Delta p}_i^{\mathrm{A}\mathrm{G}\mathrm{C}}$ AGC reference for the ith area

${\Delta p}_i^{\mathrm{W}}$ Deviation of wind power in the ith area

${\Delta p}_i^{\mathrm{W},\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}$ Predicted wind power in the ith area

${\xi }_i$ Prediction error of wind power in the ith area

$p_i^{\mathrm{E},\mathrm{d}},p_i^{\mathrm{E},\mathrm{c}}$ Discharge and charge power of ith ESS cluster

$p_k^{\mathrm{d}},p_k^{\mathrm{c}}$ Discharge and charge power of kth ESU

$s_i^{\mathrm{E}}$ Energy level of ith ESS cluster

$s_k$ Energy level of kth ESU

${SoC}_k$ SoC level of kth ESU

${\mathit{r}}_{\mathit{0},t}$ Base control output of the ESS

K Gain matrix of feedback control

Parameters

$H_i$ Equivalent system inertia

$D_i$ Equivalent system damping

$T_i^{\mathrm{t}}$ Turbine time constant

$R_i$ Droop coefficient

$T_{ij}$ Synchronizing coefficient

$B_i$ Frequency bias factor

$K_i^{\mathrm{I}}$ Integral coefficient

${\eta }_k^{\mathrm{c}},{\eta }_k^{\mathrm{d}}$ Charge and discharge efficiencies

$P_k^{\max }$ Rated power

$S_k^{\max }$ Rated energy capacity

${\Delta f}^{\max }$ Bound of frequency

$\underset{¯}{SoC},\overline{SoC}$ Lower and upper bound of SoC

${\mathit{D}}_i$ Communication matrix of ith ESS cluster