Real-Time Flexibility Assessment for Power Systems with High Wind Energy Penetration

Abstract

To reduce the reliance on fossil fuel-based generation, many countries expand the use of renewable energy sources (RES) for electricity production. The stochastic and intermittent nature of such sources (i.e., wind and solar) poses challenges to the stable and reliable operation of the electric power system (EPS) and requires sufficient operational flexibility. With continuous and random changes in the EPS operational conditions, evaluating the system flexibility in a standardized manner may improve the robustness of planning and operating procedures. Therefore, the development of fast algorithms for determining system flexibility is a critical issue. In this paper, the flexibility of the EPS with high wind energy penetration is calculated in real time. In this context, the EPS flexibility is understood as the ability of the system to maintain a balance under irregular and short-term active power variations during a specified time by using the flexibility resources. The EPS flexibility calculation relies on a deterministic method developed to qualitatively and quantitatively assess the EPS readiness to changes in load. Accurate wind power forecasts and the observance of the electric circuit law when solving the optimization problem allow for determining the actual value of the EPS flexibility during a considered time.

1. Introduction

With the government targets set to reduce the dependency on fossil fuel energy sources, there is a rapid increase in using renewable energy sources (RES) for electricity generation. The most abundant renewable energy sources, wind and solar, are well received for deployment in many existing power systems (for example, Denmark, Ireland, and Germany) [1]. These sources are stochastic and intermittent in nature, which implies that in addition to traditionally variable demand, the power systems have to deal with rising variabilities in the power supply [2].

For the efficient management of electric power systems, including those with RES, forecasting of electric demand and wind and solar power outputs is utilized. Based on such forecasts, scheduling of controllable and flexible sources of energy (e.g., conventional power plants, batteries, etc.) is performed. The actual operation of the system elements, however, differs from the predicted (planned) one due to forecast errors, which consist of errors in the forecasting method, an unexpected change in electric demand, and an unexpected change in the generation at a wind (solar) power output. Forecast errors cannot be predicted, therefore, they belong to the category of uncertainties in the form of irregular and short-term changes in active power. While increasing RES penetration, it is essential to take measures to improve the system flexibility and prevent future stability and reliability issues.

The stable and safe operation of the electric power system (EPS) is possible if there are reserves to maintain the power balance between the generation and consumption under normal and post-emergency conditions. The contingency reserves utilized under normal operating conditions (to perform regulation services) can pose a threat to the system functioning. This situation, however, can be avoided. To this end, the imbalances caused by uncertainties should be compensated for at the expense of the flexibility sources, i.e., the power regulation reserves that are not planned to be used under emergencies. The flexibility sources include [3] the power system reserve [4,5], demand-side management [6,7], energy storage systems, dispatchable distributed generators, grid interconnection [8], and multi-energy systems [9].

The problems associated with the flexibility of power systems are solved by researchers in many countries. They have developed probabilistic and deterministic methods for calculating flexibility, as well as methods based on artificial intelligence.

A unified four-element framework (time, uncertainty, action, and cost) for determining the power system flexibility is presented in [10]. Based on this framework, the authors of [10] have developed a flexibility metric that factors in the system operation and transmission network constraints. The proposed metric quantifies the broadest variation range of uncertainty the power system can meet. The robust optimization technique is used to compute this metric [10].

The study in [11] examines the power system flexibility that becomes available via demand response. The authors consider the thermostatically controlled loads (TCLs) as an important class of demand response. The proper coordination of TCLs plays a significant role in providing different services to the grid (renewable integration, peak shaving, etc.). The control model of TCLs is integrated into the system-level operation and control is calculated. The set of valid power profiles of individual TCLs is represented by a polyhedron. The Minkowski sum is calculated to determine the aggregated flexibility. The study [11] relied on an optimization algorithm for approximating a polynomial by homotheties of a given convex set represented by a virtual battery model.

The insufficient ramping resource expectation (IRRE) metric is developed to measure the flexibility of a power system [12]. It is used to recognize the time intervals during which the power system experiences a shortage of flexible resources, and, consequently, there is a probability that it will not be able to cope with the predicted or unpredicted changes in the net load. The value of the abovementioned probability is calculated from the cumulative probability. The authors [12] give a detailed description of the algorithm for calculating IRRE for the time horizons of interest.

In [13], the flexibility metric is developed and used to accurately estimate the flexibility level of a power system for various time horizons and directions. This metric is based on a kernel density estimator and means that the EPS has insufficient ramp resources if the probability of the flexibility residual is less than zero. The flexibility residual is considered to be a mismatch between the available flexibility and the forthcoming net load ramps. The kernel density estimator needs to estimate the probability density function of the time series, where the Gaussian kernel is chosen as the kernel function. The authors of [14] simulate multiple scenarios considering different demand and generation profiles to calculate a flexibility metric.

The papers [15,16,17,18,19,20] are devoted to the use of an artificial neural network (ANN) to calculate the flexibility of an EPS. In [15,16,17], an ANN is used to model the flexibility of distributed energy sources. The trained ANN assesses whether the flexibility of specific distributed energy resources is available for a given load curve. The advantage of this approach is that there is no need to create accurate models of the allocated resources. The authors [18] explore and define the flexibility of distributed energy resources. In [19] a new approach for quantifying the network flexibility potential is introduced. In that work, historical demand data are classified using an ANN.

Many of the abovementioned methods share some similarities where in order to calculate the EPS flexibility, the simulation of multiple scenarios has to be performed.

In this research, the flexibility of an EPS with a wind farm (i.e., the ability of the system to maintain a balance with irregular and short-term changes in active power during a specified time at the expense of flexibility sources) is calculated in real time for n minutes ahead. Since the flexibility metric is assumed to be calculated in real time using the current state variables and very short-term forecast variables, there is no need to simulate any scenario. The focus is on forecasting the power generated by a wind farm. There is a whole host of studies on the creation of an ANN for wind power forecasting in the world [21].

This paper shows the significance of assessing the flexibility of an electric power system with a wind farm and a battery in real time for n minutes. The battery is one of the essential components for providing the flexibility of such power systems. The EPS flexibility calculation involves the development of a deterministic method to qualitatively and quantitatively assess the EPS readiness to the load variation. A criterion is proposed for determining the combination of imbalances arising from forecast errors (uncertainties), which, when exceeded, make the EPS incapable of maintaining the power balance (cause flexibility losses). This problem is solved as an optimization problem with equality and inequality constraints. To check if the constraints, which are the electrical circuit equations, are met, a load flow solution is found. The optimization parameters are the components of vector z, which are used for modeling a variety of loads.

The main contributions of this paper are as follows:

• A methodology for real-time flexibility assessment for a power system with high wind energy penetration is developed. This methodology uses measurements as input data, handles them, and then provides information about the EPS flexibility for n minutes ahead, taking into account the battery control strategy;
• Wind power forecasting integrated within the methodology is performed using state of the art forecasting method (variational mode decomposition and long short-term memory (VMD-LSTM));
• A new term “operating state with minimum flexibility” is introduced.

The paper has the following structure. The second section presents the purpose and concept of the research. The third section focuses on the EPS flexibility calculation and describes the modeling of the EPS flexibility sources. This section also presents, in a general form, information related to the processing of measurement data (data collection system and state estimation); and makes a forecast of the active power at the wind farm. The same section gives the main idea of the method, presents the optimality criterion and the objective function for calculating the maximum loads, and sets the constraints. The fourth section discusses the research results. The conclusion contains the analysis of the work performed.

2. Purpose and Concept of the Research

The flexibility assessment is one of the most significant issues in the operational dispatch control of an EPS with renewable energy sources. Current information on the EPS flexibility can become critical when making a decision on adjusting the operating states in unforeseen situations of the EPS functioning.

This research proposes calculating the flexibility with a method based on the statement that the EPS has flexibility if the power balance is maintained during the considered time under uncertainties on the part of power consumers and producers. First, we determine the combinations of loads for the given values of generating capacities (dispatch solution), in which the balance is ensured by putting into operation the entire power reserve intended to smooth out the uncertainties. This EPS operating state is characterized by the minimum flexibility (hereinafter referred to as the operating state with minimum flexibility) since an insignificant excess of load at a node will lead to a power shortage in the EPS or to the use of an emergency reserve. The latter is unacceptable as it will pose a threat to the EPS’s safe operation. Then the flexibility metric is calculated.

The power system flexibility in a normal operating state at the time t + n moment, with known information about the configuration and parameters of the EPS equivalent circuit (in Figure 1. EPS description), is calculated according to the following methodology (Figure 1):

• Collect information about the operating parameters, perform a state estimation (blocks I, II);
• Specify the scheduled operating state according to the measurement information (blocks III, IV). To do this, forecast the power at the wind farm using an ANN and find a load flow solution;
• Calculate the flexibility of the refined planned operating state (blocks V, VI). First, determine the operating state with minimum flexibility. This means calculating the steady-state variables (SSV) that include the largest loads for a given generation. This is an optimization problem with equality and inequality constraints.

Further, the process of operating a state flexibility assessment at the considered time (t: t + n) is described in detail using the developed algorithm.

3. Real-Time Flexibility Assessment

3.1. Data Acquisition

Traditionally, electricity is generated by large power plants and delivered to the end user via transmission and distribution networks [22]. These networks are equipped with measurement devices (e.g., remote terminal units (RTUs), intelligent electronic devices (IEDs), and phasor measurement units (PMUs)) connected to a local area network (LAN) along with supervisory control and data acquisition (SCADA) systems [23].

There are several sources of measurements in the EPS (Figure 2):

• SCADA systems, IEDs, PMUs, and RTUs provide real-time voltage, current, and power flows measurements;
• Meter data management systems (MDMS) along with a customer information system (CIS) are used to monitor the network consumption and distributed energy resources (DER) output. In case of insufficient data redundancy, a pseudo measurement generator is used;
• Automated mapping and facility management (AM/FM) systems along with geographic information and outage management systems (GIS and OMS, respectively) monitor the equipment connectivity status.

3.2. State Estimation

The state estimator is the data processing algorithm used to filter raw measurements from noise, detect and eliminate gross errors, and determine the optimal estimates for the system state. It includes topology processing, an observability analysis, a state estimation, and bad data detection [24,25].

The measurement set used by the state estimator looks as follows:

$$\overline{y} = \left(U_i,P_i,Q_i,P_{ij},Q_{ij},{\delta }_i\right),$$

(1)

where $P_i = P_{l\left(i\right)}+P_{g\left(i\right)}$, $Q_i=Q_{l\left(i\right)}+Q_{g\left(i\right)}$, $U_i$ is the magnitudes of nodal voltages; $P_{g\left(i\right)}, Q_{g\left(i\right)}$ are the active and reactive power generation at the nodes; $P_{l\left(i\right)}, Q_{l\left(i\right)}$ are the active and reactive power load at the nodes; $P_{ij}, Q_{ij}$ are power flows in the transformers and the lines; ${\delta }_i$ is the voltage phases at the nodes.

The objective function used in the traditional state estimation is as follows:

$$J\left(x\right) = {(\overline{y}-y\left(x\right))}^T{R}_y^{-1}\left(\overline{y}-y\left(x\right)\right),$$

(2)

where $R_y$ is a covariance matrix of measurement errors (diagonal matrix), $x = \left(U,\delta \right)$ is the state vector (a part of the state variables, which is used to find the load flow solution).

The state vector calculation involves the minimization of criterion Equation (2) with respect to the voltage phase and magnitude.

$$J\left(x\right)\to min.$$

(3)

The derivative $\partial J\left(x\right)/\partial x$ is set to zero, and the system of non-linear equations is solved:

$$H^T{R}_y^{-1}\left(\overline{y}-y\left(x\right)\right) = 0$$

(4)

subject to

$$w\left(y,y1\right) = 0,$$

where $H = \partial y/\partial x$ is a Jacobi matrix, $w\left(y,y1\right) = 0$ represents the equations of the electric circuit connecting the measured (y) and unmeasured (y1) variables.

The system Equation (4) is linearized at each iteration i. The system of linear equations is solved, where the correction vector is calculated by the equation:

$$\Delta x^{\left(i\right)} = {\left(H^{T\left(i\right)}{R}_y^{-1}{H}^{\left(i\right)}\right)}^{-1}(H^{T\left(i\right)}{R}_y^{-1}\left(\overline{y}-y\left(x^{\left(i\right)}\right)\right).$$

(5)

The components of the state vector are calculated by the equation:

$$x^{\left(i+1\right)}=x^{\left(i\right)}+\Delta x^{\left(i\right)}.$$

(6)

If $\Delta x<{\epsilon }_x$, the iteration process converges, where ${\epsilon }_x$ is the accuracy of the iterative process. Then, a load flow solution is obtained using the state vector.

3.3. Wind Power Forecasting Using VMD-LSTM

The wind farms power outputs have both fluctuating and random components. Therefore, forecasting of the future wind energy production might be a challenging task. To adequately determine the EPS flexibility, it is particularly important for system operators to have as accurate forecasts of power outputs as possible.

The forecasting technique used in this paper is based on the hybrid variational mode decomposition and long short-term memory (VMD-LSTM) artificial neural network [26]. Hybrid models can recognize hidden signal features in time series (including wind power time series) and learn from them. The original wind power time-series signal is decomposed into multiple sub-series, which are less chaotic and, consequently, more predictable. For this purpose, variational mode decomposition (further VMD method) is used. The VMD method is described in detail in [27]. To implement this method, the following mathematical functions are used: Wiener filtering, Hilbert transform and frequency mixing, and heterodyne demodulation.

The wind power time series are decomposed according to the physical signal properties. With such a decomposition, the original wind power signal is divided into three components: the long-term trend, short-term fluctuation, and ultra-short randomness.

The algorithm for forecasting wind power consists of five steps, described below.

• Select an ANN type. This study uses an LSTM ANN. It is a kind of recurrent neural network that overcomes the memory limitations by adding particular cells aimed at filtering new information (preserve or ignore);
• Form a data set for training, validation, and testing. In this study, historical wind power measurements are used as a data set;
• Initialize LSTM. Determine the LSTM starting parameters. The starting parameters are as follows: ${\sigma }_s$ and ${\sigma }_t$ are activation functions (${\sigma }_s$ denotes the state (sigmoid), and ${\sigma }_t$ denotes the gate (hyperbolic tangent)); W, M, and b are input weights, the recurrent weights, and the bias of each component, respectively; and $c_t$ is the cell memory state, $h_t$ is the hidden state, $x_t$ is a time-series vector.
• Train LSTM. Use the method of backpropagation for training LSTM. The key element of LSTM is a cell state. During the LSTM training, the memory cell state changes. LSTM is considered to be trained if the memory cell stops changing. The following equations describe the forward pass of an LSTM unit with the input, output, and forget gates and with the memory cells at each time step t (Figure 3):

The sigmoid layer of the LSTM cell determines what information should be removed from the cell of state

$$f_t={\sigma }_s\left(W_f{x}_t+M_f{h}_{t-1}+b_f\right).$$

(7)

The input gate determines what part of the new information needs to be added to the memory

$$i_t = {\sigma }_s\left(W_i{x}_t+M_i{h}_{t-1}+b_i\right).$$

(8)

The tanh layer creates the cell candidate

$$g_t = {\sigma }_t\left(W_g{x}_t+M_g{h}_{t-1}+b_g\right).$$

(9)

The new information is added

$$c_t = f_t{c}_{t-1}+i_t{g}_t$$

(10)

The output gate defines the extent to which the existing part of the memory contributes to the output

$$o_t = {\sigma }_s\left(W_o{x}_t+M_o{h}_{t-1}+b_o\right).$$

(11)

The hidden state is:

$$h_t = o_t{\sigma }_t\left(c_t\right).$$

(12)

5.

Test LSTM.

Figure 3. Long short-term memory.

Figure 3. Long short-term memory.

Figure 4 shows a schematic representation of the VMD-based LSTM wind power forecasting where two measurements are used to predict one point.

The wind power forecasting process may be divided into several steps.

Step 1. Select the input data. Use several previous measurements (two in Figure 4) to predict the wind power output t + 1 min ahead.

Step 2. Perform the VMD to decompose each input datum into three features. These features are as follows: the long-term trend ($P_{t-1}^{rand}$, $P_t^{rand}$), short-term fluctuation ($P_{t-1}^{fluc}$, $P_t^{fluc}$), and ultra-short randomness ($P_{t-1}^{trand}$, $P_t^{trand}$).

Step 3. Forecast wind power. Three components of wind power are predicted using an LSTM artificial neural network.

Step 4. Build the predicted point $P_{t+1}$ from $P_{t+1}^{rand}$, $P_{t+1}^{fluc}$, and $P_{t+1}^{trand}$.

3.4. Computing Load Flow Solution t + n Minute Ahead

The load flow calculation aims to clarify the scheduled state variables at time t + n according to the information available at time t. In this case, the scheduled state variables are referred to as the load patterns and generation profiles (unit commitment solutions), to which the time instant t + n belongs. In the load flow calculation, injections of active and reactive power at all the nodes, except for the slack one, are used as independent state variables. The variables are refined in the following order:

• The wind farm injections are predicted for n minutes ahead;
• The capacity values for the load and generator units are taken from the data of the scheduled operating state;
• The battery power is calculated according to the control strategy for the battery energy storage system (BESS) based on the forecast data on the wind farm power and the scheduled load data. The BESS control strategy can vary. For example, the strategy can suggest the use of a BESS to maximize the renewable energy penetration or to smooth a profile of power output from the wind and solar power plants. This study uses the former [28];
• The reactive power values are calculated based on the rule of maintaining the proportion between active and reactive power;
• The initial voltage approximations are set. The initial approximations are the voltages calculated with the state estimation at time t;
• Based on the known active and reactive injections at the nodes, the load flow solution is computed.

3.5. The Flexibility Assessment

The power system flexibility is characterized by the following three indices: ramping limit (ΔR), power capacity (P), and energy capacity (E). Each index is a time derivative of the previous index [3].

The calculation of the flexibility of the EPS resources:

• Assumes that the ramping limit of the flexible source (ΔR) is equal to the lead time for which the flexibility is calculated;
• Factors in the maximum and minimum generation (P) index;
• Takes into account that the final energy capacity of the battery (E) limits its capabilities.

3.5.1. Modeling the Flexibility of EPS Components

Model of generator flexibility of a conventional power plant

The flexibility available from each generator $F_g$ is determined by the power that can be additionally generated over the time under consideration and is calculated using the equation:

$$F_g=P_{g\left(\mathrm{per}\right)}-P_{g\left(c\right)},$$

(13)

where $P_{g\left(\mathrm{per}\right)}$ is the regulating reserve (employed in normal condition), $P_{g\left(c\right)}$ is the active power of the generator at the moment in question.

Model of battery flexibility

The flexibility available from the battery $F_B$ is determined by the state of charge and the current power. If the battery is charged within the specified limits:

$$So{C}_{min}<SoC\left(t\right)<So{C}_{max},$$

(14)

then the power output is calculated using the equation:

$$F_B=P_{B\left(max\right)}-P_{B\left(c\right)},$$

(15)

otherwise,

$$F_B=0,$$

(16)

where SoC is the state of charge, $So{C}_{min}$ is the minimum level of charge, $So{C}_{max}$ is the maximum level of charge, $P_{B\left(c\right)}$ is the active power of the battery at the moment in question.

Model of system flexibility at the time in question

The power system flexibility at the time in question ($F_{Sys}$) is determined as a sum of flexibilities of all the facilities in the system. When two kinds of flexibility are available in the system, $F_{Sys}$ is calculated by the equation:

$$F_{Sys}={{\displaystyle \sum }}_1^{n1}{F}_g+{{\displaystyle \sum }}_1^{n2}{F}_B,$$

(17)

where n1 is the number of generators at conventional plants, n2 is the number of batteries.

3.5.2. Metric of Flexibility

Whether the power system is flexible or not during the time [t: t + n] is determined by calculating a flexibility metric:

$$F_S=P_{calc}^{max}-P_L,$$

(18)

where $P_{calc}^{max}$ is the largest active loads when the power system can remain flexible, $P_L$ is the active powers in a scheduled operating state.

$$P_{calc}^{max}={{\displaystyle \sum }}_{i=1}^{n3}{P}_i\left(z_i\right),$$

(19)

$$P_L={{\displaystyle \sum }}_{i=1}^{n3}{P}_{L\left(i\right)},$$

(20)

where $P_i\left(z_i\right)$ is the maximum load at node i, at a set loading of generators (at a given dispatch solution), $n_3$ is the number of load nodes.

If F_S > 0, then the considered EPS is flexible (during the time [t: t + n]). The value $F_S$ is the largest range of uncertainty, within which there is no danger for the power system. $F_S$ shows how much flexibility the power system has.

3.5.3. Determination of Operating State with Minimum Flexibility. Objective Function

The operating state with the minimum flexibility is calculated as a problem of nonlinear optimization (we have non-linear constraints). The optimality criterion is the maximum difference between the values of active powers at nodes with uncertainty in the scheduled operating state and in the operating state with minimum flexibility. The optimization parameters are the components of vector z, which are responsible for modeling the load value at load nodes from the minimum to maximum values.

The objective function, i.e., the maximum sum of differences between the modeled and scheduled powers at nodes with uncertainty, is written as follows:

$$J\left(z\right)={{\displaystyle \sum }}_{i=1}^r\left(P_i\left(z_i\right)-{\overline{P}}_i\right)={{\displaystyle \sum }}_{i=1}^r\Delta P_i\left(z_i\right)\to max,$$

(21)

where r is the number of nodes with uncertainty. In this case, these are the load nodes, nodes where control actions are performed, and a slack node, $P_i\left(z_i\right)$ is calculated as follows [10]:

$$P_i\left(z_i\right)=z_i{P}_i^{min}+\left(1-z_i\right)P_i^{max}.$$

(22)

A distinctive feature of the research is the method of forecasting the active power of a wind farm, which gives more accurate results in comparison with load forecasts. Therefore, wind farm power is considered to be an exact value and is not entered into the objective function. The predicted battery power value calculated based on the predetermined battery control strategy is also not input to the objective function. The uncertainties associated with an error in a forecast of power at the wind farm and battery are compensated for at the balancing power plant.

Each component Equation (21), given Equation (22), can be represented as follows:

$$J\left(z\right)=\Delta P_i\left(z_i\right)=\left(P_i\left(z_i\right)-{\overline{P}}_i\right)=P_i^{max}-z_i{P}_i^{max}+z_i{P}_i^{min}-{\overline{P}}_i.$$

(23)

When all the values independent of the optimization parameters are excluded, the objective function has the form:

$$\underset{z}{\max }{{\displaystyle \sum }}_{i=1}^r-(P_i^{max}-P_i^{min})z_i \forall z\in \left[0, 1\right].$$

(24)

In Equations (23) and (24), i is the number of a node with uncertainty, ${\overline{P}}_i$ is the value of active power at node i in the scheduled operating state; $P_i\left(z_i\right)$ is the relationship between active power and value z responsible for the variation in values of power at node i, $P_i^{max}$, $P_i^{min}$ are the maximum (minimum) values of consumed (generated) power at node i. The maximum generated power is the maximum power that can be used in the normal operating state, in particular, to compensate for the uncertainties.

The search for optimal parameters is performed given the equality and inequality constraints. The equality constraints are the balance of active power in system Equation (25) and the equation that controls the implementation of the electric circuit laws Equation (26). The inequality constraints are introduced to control the transfer capability of the lines and the power limits at nodes with uncertainty Equations (27)–(29).

The constraints have the following form:

$$\Delta P_s={{\displaystyle \sum }}_{j=1}^{n-r}{\overline{P}}_j+{{\displaystyle \sum }}_{i=1}^r{P}_i\left(z_i\right)+\pi =0,$$

(25)

$$P_i\left(z_i\right)-P_i\left(U_i,{\delta }_i\right)=0,$$

(26)

$$P_{l-j}<P_{l-j}^{max},$$

(27)

$$P_i^{min\left(CA\right)}<P_i^{CA}<P_i^{max\left(CA\right)},$$

(28)

$$P_i^{min\left(S\right)}<P_i^S<P_i^{max\left(S\right)},$$

(29)

where n is the number of nodes in the grid; $\Delta P_s \mathrm{is} \mathrm{the} \mathrm{active} \mathrm{power} \mathrm{balance} \mathrm{in} \mathrm{the} \mathrm{system}; P_i\left(z_i\right)$ is the node i power value obtained from solving the optimization problem; $P_i\left(U_i,{\delta }_i\right)$ are state variables at node i; $U_i,{\delta }_i$ are the voltage magnitude and phase at node i; $P_{l-j}^{\mathit{max}}$ is the transfer capability limit of the transmission line limited by nodes l, j; $P_i^{min\left(CA\right)}$, $P_i^{max\left(CA\right)}$, $P_i^{min\left(S\right)}$, $P_i^{max\left(S\right)}$ are the maximum and minimum active power at the controlled and slack nodes, respectively.

The power balance in system Equation (25), given Equation (22) and after some transformations, is written as follows:

$${{\displaystyle \sum }}_1^r{B}_i{z}_{i }={{\displaystyle \sum }}_1^r{P}_i^{max}+{{\displaystyle \sum }}_1^{n-r}{\overline{P}}_j+\pi ,$$

(30)

where n is the number of nodes in the system, $\pi$ is the system active power losses calculated in the operating state. $B_i=(P_i^{max}-P_i^{min})$ is for the load nodes, $B_i=(-P_i^{max}+P_i^{min})$ is for the generator nodes.

Equation (26) guarantees the observance of Ohm’s and Kirchhoff’s laws since it includes the steady-state variables. To calculate the steady-state variables, a system of nodal equations is built in the form of power balances:

$$U_l^2{G}_l+{{\displaystyle \sum }}_{j=1}^m{P}_{l-j}\le {\xi }_P,$$

(31)

$$U_l^2{D}_l+{{\displaystyle \sum }}_{j=1}^m{Q}_{l-j}\le {\xi }_Q,$$

(32)

where $G_l$, $D_l$ are the active and reactive conductivity of the shunt at node l; $P_{l-j}$, $Q_{l-j}$ are the active and reactive power in the lines limited by nodes l, j; m is the number of adjacent nodes. This system of equations is solved by Newton’s iterative method, where ${\xi }_P$, ${\xi }_Q$ is the accuracy of the iteration process. The achieved accuracy of Equations (31) and (32) means the end of the iterative process.

4. Case Study

4.1. Description of a Scheme and a Scenario

The developed method is applied to calculate the power system flexibility in the normal operating state and is verified using the test scheme shown in Figure 5:

The scheme consists of six tie lines and six nodes, including a wind farm (node 1), a battery (node 5), conventional plants (nodes 2 and 6), and loads (nodes 3 and 4). Nodes 2, 3, 4, and 6 are assigned to be the nodes with uncertainty. Nodes 1 and 5 are considered to be the nodes without uncertainties since the uncertainties associated with the wind farm and the battery are assumed to be compensated for at the balancing power plant (node 6).

The calculations are performed according to the scenario: calculate the power system flexibility in two normal operating states 10 min ahead. The operating states differ from each other by the power of the battery. The load uncertainties are compensated for at the controlled plant. The uncertainties exceeding the permissible size of the control action are compensated for by the power at the balancing power plant.

The ramping limit (ΔR) is assumed to be 10 min. This means that in 10 min, the entire reserve available at a conventional plant, designed to compensate for uncertainties, is turned on, and the battery delivers the maximum power. The maximum and minimum values of active power at nodes with uncertainty and the active powers of two operating states are shown in

Table 1. Initial information (MW).

Number of the Node	${\mathit{P}}^{\mathit{m}\mathit{a}\mathit{x}}$	${\mathit{P}}^{\mathit{m}\mathit{i}\mathit{n}}$	Operating States
			1	2
1			16.14 (16.04)	15.12 (15.0)
2	30	21	22.98	21.52
3	−32	−26	−24.59	−23.05
4	−22	−16	−18.39	−19.63
5			0.05	6.39
6	10	6	7.8	4.67

. The powers are taken from the operating state scheduled 10 min ahead.

4.2. Flexibility Calculation for a Six-Node System

This problem-solving process can be represented by several steps.

Step 1. Acquire the data. Perform the state estimation.

Step 2. Forecast the wind power. Make a very short-term forecast (10 min ahead). Use the historical wind power measurements as a training data set. The answers of trained LSTM are 16.14 and 15.12. (The real values of wind power output are 16.04 and 15.0).

Step 3. Determine the state variables 10 min ahead. The results are shown in Table 1, columns 4 and 5. The real values of wind power output are shown in brackets.

Step 4. Describe the objective function and constraints following the requirements for the input data of the software programs used.

The objective function

$$\left(P_2\left(z_2\right)-{\overline{P}}_2\right)+\left(P_3\left(z_3\right)-{\overline{P}}_3\right)+\left(P_4\left(z_4\right)-{\overline{P}}_4\right)+\left(P_6\left(z_6\right)-{\overline{P}}_6\right)\to max,$$

(33)

is written as follows:

$$\underset{z}{\max }{{\displaystyle \sum }}_{i=1}^4(-B_i{z}_i) , \forall z\in \left[0, 1\right],$$

(34)

where $B_i=P_i^{max}-P_i^{min}$, i = 2, 3, 4, 6.

The power balance in the system is written in the following way:

$${{\displaystyle \sum }}_{i=1}^4(B_i{z}_i)=P_3^{max}+P_4^{max}-P_2^{max}-P_6^{max}-{\overline{P}}_1-{\overline{P}}_5-\pi .$$

(35)

The constraints that control the power balances at the nodes have the form

$$P_2\left(z_2\right) - P_2\left(U_2,{\delta }_2\right) \le \xi ,$$

(36)

$$P_3\left(z_3\right) - P_3\left(U_3,{\delta }_3\right) \le \xi ,$$

(37)

$$P_4\left(z_4\right) - P_4\left(U_4,{\delta }_4\right) \le \xi ,$$

(38)

$$P_6\left(z_6\right) - P_6\left(U_6,{\delta }_6\right) \le \xi ,$$

(39)

where $z_2,z_3, z_4,z_6$ are optimization parameters, $P_2\left(U_2,{\delta }_2\right), P_3\left(U_3,{\delta }_3\right)$, $P_4\left(U_4,{\delta }_4\right)$, $P_6\left(U_6,{\delta }_6\right)$ are active powers at the nodes, which are the steady-state variables. The constraints (31) and (32) are factored in in the load flow calculation.

The constraints of the form Equations (27)–(29) are built relying on the data from Table 1.

Step 5. Solve the optimization problem considering the set constraints. The Matlab environment is employed with the function fmincon. The result of the optimization problem is $z_2,z_3, z_4,z_6$. Equation (22) is then used to compute the optimal steady-state variables, which are the variables of the operating state with minimum flexibility.

Step 6. Assess the flexibility. Equation (18) is used to measure the flexibility of the power system for 10 min.

4.3. Analysis of Results and Discussion

The results are summarized in

Table 2. Calculation results (P (MW)).

Number of the Node	Operating State 1		Operating State 2
	${\mathit{P}}_{\mathit{i}}^{\mathit{c}\mathit{a}\mathit{l}\mathit{c}}\left(\mathit{z}\right)$	${\mathit{F}}_{\mathit{S}\left(\mathit{i}\right)}$	${\mathit{P}}_{\mathit{i}}^{\mathit{c}\mathit{a}\mathit{l}\mathit{c}}\left(\mathit{z}\right)$	${\mathit{F}}_{\mathit{S}\left(\mathit{i}\right)}$
1
2	29.2		29.6
3	−28.0	3.41	−32.0	8.95
4	−22.0	3.61	−22.0	2.37
5
6	9.69		9.12
$F_S$		7.02		11.32

. Columns 2 and 4 contain the active power values determined in accordance with Equation (22), where z is the result of optimization. Columns 3 and 5 show flexibility $F_{S\left(i\right)}$. The last line indicates the total EPS flexibility calculated according to Equation (18).

Figure 6 shows the maximum power values (Table 1, column 2) and the optimization results (Table 2, columns 2 and 4).

The analysis of the results presented in Table 2 and Figure 6 shows that:

• A comparison of the preset maximum loads (Table 1, italics) with the maximum possible loads (Table 2, bold type) demonstrates: 1) the calculated loads do not exceed the maximum values; 2) if we take the value of loads in the scheduled operating state as 100%, then, the maximum load deviations, given the constraints set, are 13.87% and 19.63% in operating state 1 (38.82% and 12.07% in operating state 2) at nodes 3 and 4, respectively;
• The flexibility of the considered system is equal to 7.02 MW and 11.32 MW for operating state 1 and operating state 2, respectively. This means that when the battery delivers the maximum power, the highest load deviation, provided that the balance can be restored, is 11.32 MW; when the battery supplies the minimum power, the maximum load deviation is 7.02 MW.

The flexibility was calculated for 25 snapshots (sets of measurements taken simultaneously). In this case, the measurements were taken 25 times every 10 min. Figure 7 shows the active power at all the nodes in 25 normal operating states. Figure 8 shows the active power at all the nodes in the operating states with minimum flexibility. Figure 9 shows the difference between two operating states, where the difference between the active power at nodes 3 and 4 is used to calculate the flexibility ($F_s$). The analysis of the results, presented in Figure 7, Figure 8 and Figure 9, shows how the EPS flexibility differs at different time.

The developed methodology provides information about the upward flexibility. In some cases, it is important for system operators to have the knowledge of the downward flexibility, which can be a subject of further research.

In this work, the utilized BESS control strategy was used for illustrative purposes. To determine the effect of implementing different control strategies on the EPS flexibility, further research is needed.

5. Conclusions

The paper presents a method for the real-time calculation of the flexibility of an electric power system with a wind farm. An optimality criterion has been developed to identify the combination of loads that are the maximum for the planned generating capacities in the period at issue. The problem of flexibility calculation is solved as an optimization problem with a linear objective function and equality and inequality constraints. The developed method enables a quantitative metric of flexibility to be calculated.

When solving the optimization problem, the observance of Ohm’s and Kirchhoff’s laws is guaranteed by introducing a constraint as a function representing the load flow solution. The calculation of the steady-state variables involves the construction of a system of nodal equations, which is solved by Newton’s method. The state variables are maintained within the specified limits by introducing inequality constraints.

The size of flexibility is calculated for the normal operating state disregarding emergencies in the EPS. At the same time, when using the obtained information to adjust the state variables, the system is guaranteed not to be exposed to the threat of a system failure during the considered time. The active power balance is maintained at the expense of the flexibility sources. The emergency reserve remains unchanged, which is achieved by correctly set power constraints at the slack node (the available power is reduced by the emergency reserve magnitude).

The wind farm active power has been forecast. The developed method for its forecasting is sufficiently accurate (the error of 0.6% in the first operating state, and 0.8% in second (Table 1)). The ability of the forecasting system to predict wind power has been improved by decomposing the original wind power measurements into three components: the long-term trend, fluctuation component, and signal randomness.

The flexibility of a six-node system, which includes a wind power plant, a storage battery, two loads, and two generating plants, has been calculated. A comparative analysis of the flexibility values calculated for the different operating states shows that the EPS flexibility value depends on the state variables at the considered time instant. This confirms the need to determine flexibility in real time for some time in advance.