Optimized Auxiliary Frequency Control of Wind Farm Based on Piecewise Reduced-order Frequency Response Model

Xu Zhang; Chen Zhao; Junchao Ma; Long Zhang; Dan Sun; Chenxu Wang; Yan Peng; Heng Nian

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

Optimized Auxiliary Frequency Control of Wind Farm Based on Piecewise Reduced-order Frequency Response Model PDF

- ORCID：
Xu Zhang ¹
✉
- ORCID：
Chen Zhao ¹
✉
- ORCID：
Junchao Ma ²
✉
- ORCID：
Long Zhang ¹
- ORCID：
Dan Sun ¹
✉
- ORCID：
Chenxu Wang ²
✉
- ORCID：
Yan Peng ²
✉
- ORCID：
Heng Nian ¹
✉

1. School of Electrical Engineering, Zhejiang University, Hangzhou 310058, China； 2. Electric Power Research Institute of State Grid Zhejiang Electric Power Corporation, Hangzhou 310014 China

Updated：2024-05-20

DOI：10.35833/MPCE.2023.000448

OUTLINE

Abstract

With the increasing wind power penetration in the power system, the auxiliary frequency control (AFC) of wind farm (WF) has been widely used. The traditional system frequency response (SFR) model is not suitable for the wind power generation system due to its poor accuracy and applicability. In this paper, a piecewise reduced-order frequency response (P-ROFR) model is proposed, and an optimized auxiliary frequency control (O-AFC) scheme of WF based on the P-ROFR model is proposed. Firstly, a full-order frequency response model considering the change in operating point of wind turbine is established to improve the applicability. In order to simplify the full-order model, a P-ROFR model with second-order structure and high accuracy at each frequency response stage is proposed. Based on the proposed P-ROFR model, the relationship between the frequency response indexes and the auxiliary frequency controller coefficients is expressed explicitly. Then, an O-AFC scheme with the derived explicit expression as the optimization objective is proposed in order to improve the frequency support capability on the premise of ensuring the full release of the rotor kinetic energy and the full use of the effect of time delay on frequency regulation. Finally, the effectiveness of the proposed P-ROFR model and the performance of the proposed O-AFC scheme are verified by simulation studies.

Keywords

Frequency response; piecewise reduced-order; wind farm; auxiliary frequency control

I. Introduction

WITH the deterioration of the global environment, renewable energy power generation such as wind power and photovoltaic has been paid much attention [

1]-[3]. Nowadays, more and more wind power generators with conventional control strategy are connected to the power grid through power electronic devices, which do not have the same frequency regulation capability as the synchronous generator (SG). Therefore, in order to improve the frequency stability of the power system, auxiliary frequency control (AFC) is adopted by wind farm (WF) to simulate the active support capability of the SG [4]-[6]. However, the effectiveness of the AFC depends on the design of its control coefficients [7]-[9]. Therefore, in order to maximize the frequency support capability of WF, it is necessary to reveal the relationship between the frequency response indexes and the control coefficients, and then study the optimization scheme for the AFC.

The accurate evaluation of frequency response indexes depends on the precise frequency response modeling. Existing studies on the frequency response modeling are mostly based on the single-machine equivalent model due to the merit of simple structure [

10], [11]. The single-machine equivalent model usually includes the average system frequency (ASF) model and the system frequency response (SFR) model. The ASF model has satisfactory accuracy, but the order of the model is fairly high, so it is difficult to derive the analytical expression of frequency nadir (FN) during frequency faults [12]. Based on the ASF model, the SFR model with a simpler structure is further proposed in [13]. For the second-order SFR model, the analytical expression of frequency response indexes can be derived by inverse Laplace transform. However, the accuracy of the SFR model is poor due to excessive simplification and approximation. Besides, the SFR model in [13] only takes SGs into account, which is not applicable for multi-machine power system integrated with wind power generation. In general, it is necessary to simplify the frequency response model of the multimachine system while preserving the high accuracy.

In order to reduce the order of ASF model or improve the accuracy and applicability of SFR model, some improved frequency response models have been proposed. In [

14], an improved SFR model under the joint action of SG and WF is established. However, the accuracy of the improved model is poor due to ignoring the influence of the change in the operating point of wind turbine (WT) on the output active power. In order to improve the accuracy of the SFR model, an improved SFR model considering the change in operating point of WT is further proposed in [15], but it is not applicable for the power system with multiple heterogeneous devices. Moreover, it is difficult to derive the analytical expression of the frequency response indexes such as FN in this model due to its high-order structure. In order to reduce the order of ASF model, an improved ASF model with low-order structure is achieved in [16] by weighting the parameters and ignoring the time constant with little effect on the frequency response. However, the reduced-order method of the improved ASF model is similar to that of the classical SFR model, and the over-simplification of the model reduces its accuracy to a certain extent. In order to reduce the order of the frequency response model while ensuring the modeling accuracy, a generic SFR model is proposed in [17]. However, the power disturbance experiment must be carried out in advance to obtain the frequency response data for the evaluation of model parameters, which might cause the unit off-grid or other hazards under unstable power grid. Therefore, there is still a gap on the frequency response model with all merits of high accuracy, ability to derive the analytical expression of the frequency index, and strong applicability.

After establishing the frequency response model, the explicit expression of the frequency response indexes can be derived as a function of control coefficients in WF, thus the optimization scheme of the auxiliary frequency controller can be further studied. In order to improve the frequency support capability, the maximum rate of change of frequency (RoCoF) RoCoF_max is regarded as an index for adjusting the coefficients of auxiliary frequency controller in [

18] and [19]. These schemes ignore the time delay of the active power output by WF. However, the AFC of WF includes frequency detection, communication, and instruction generation, so a time delay between frequency drop and output power response is unavoidable. Therefore, in order to accurately reflect the frequency regulation capability of WF, it is more reasonable to regard the average RoCoF RoCoF_avg over a period as the index. Moreover, the effect of time delay on frequency regulation is nonlinear, and appropriate time delay can effectively improve the frequency response characteristics. Therefore, it is necessary to consider the effect of time delay when optimizing the auxiliary frequency controller. In order to maximize the frequency support capability, some optimized frequency regulation schemes of WT utilizing rotor kinetic energy are proposed in [20] and [21]. However, the rotor speed constraint is not considered, which may cause the rotor speed to be lower than the safety value and result in a secondary frequency drop accident. In [22] and [23], the rotor speed limit is regarded as an optimization constraint to improve the frequency regulation effect of WF. However, the coupling relationship between rotor speed and system frequency is not analyzed and the quantitative formula of rotor speed is not discussed. Therefore, in order to fully release the rotor kinetic energy within the safe range, it is necessary to derive the analytical expression of the WT rotor speed.

This paper proposes a piecewise reduced-order frequency response (P-ROFR) model and an optimized auxiliary frequency control (O-AFC) scheme of the WF. The main contributions are as follows:

1) A P-ROFR model with second-order structure and high accuracy at each frequency response stage is proposed in order to simplify the full-order model. In order to improve the accuracy of the ROFR model at different frequency response stages, a P-ROFR model composed of three reduced-order models is further proposed, whose coefficients are obtained by minimizing the amplitude of the error model.

2) An O-AFC scheme with inner- and outer-layer optimizations is proposed to maximize the frequency support capability of WF. The optimization of control coefficients in the inner layer is to realize the full release of rotor kinetic energy and safe operation of the system. The optimization of time delay in the outer layer is to make full use of the effect of time delay on frequency regulation.

This paper is structured as follows. Section II introduces the full-order frequency response model for multi-machine power system. Section III proposes the ROFR model to simplify the full-order model, and further proposes the P-ROFR model to improve the accuracy at different frequency response stages. Section IV proposes the O-AFC scheme with inner- and outer-layer optimization. Section V verifies the effectiveness of the proposed P-ROFR model and the performance of the proposed O-AFC scheme by simulation studies. Finally, Section VI concludes this paper.

II. Full-order Frequency Response Model for Multi-machine Power System

When a frequency fault occurs, SGs and WFs inject active current into the power grid to support grid frequency recovery. The frequency response model for multi-machine power system can be expressed as:

Δ f = \frac{1}{2 H s + D} (- Δ P_{d} + Δ P_{W e} + \sum_{k = 1}^{N_{T}} Δ P_{G k}^{T} + {\sum_{k = 1}^{N_{H}} Δ P}_{G k}^{H})

(1)

where $Δ f$ is the system frequency deviation; $H$ is the inertia constant of the system; D is the damping factor of the system; $Δ P_{d}$ is the power disturbance; $Δ P_{W e}$ is the active power deviation output by the WF; $N_{T}$ and $N_{H}$ are the numbers of SGs in the thermal power plants and hydropower plants, respectively; and $Δ P_{G k}^{T}$ and $Δ P_{G k}^{H}$ are the mechanical power deviations output by the k^th SG in the thermal power plants and hydropower plants, respectively.

The relationship between the output power of the SG and the system frequency deviation can be obtained from the ASF model, so this section only focuses on the mathematical correlation between the output active power of the WF and system frequency. The model of the WT with AFC is shown in Fig. 1, where $P_{W e 0}$ is the initial value of the active power output by WF.

Fig. 1 Model of WT with AFC.

It can be observed from Fig. 1 that when the system frequency drops, an additional power signal $Δ P_{a d d}$ from auxiliary frequency controller in WT is added to the power signal $P_{M P P T}$ , which is output by the maximum power point tracking (MPPT) controller. According to Fig. 1, the mechanical power captured by the WT $P_{W m}$ and the active power output by the WT $P_{W e}$ can be expressed as (2). It should be noted that the AFC can be applied not only to WTs [

19], but also to WFs [23]. To reduce computational burden for steady-state and dynamic analysis, the difference in operating state among WTs is ignored, and a single-machine aggregation model for the WF is adopted in this paper [21].

\{\begin{array}{l} P_{W m} = \frac{1}{2} ρ S C_{p} V^{3} \\ P_{W e} = P_{M P P T} + Δ P_{a d d} = k_{o p t} ω_{r}^{3} - \frac{K_{d} s f + K_{p} Δ f}{T_{d} s + 1} \\ P_{W m} - P_{W e} = 2 H_{W T} ω_{r} s ω_{r} \end{array}

(2)

where $ρ$ is the air density; $S$ is the area swept by the WT; $C_{p}$ is the power coefficient; $V$ is the wind speed; $k_{o p t}$ is the MPPT coefficient; $ω_{r}$ is the rotor speed of the WT; $K_{d}$ and $K_{p}$ are the differential coefficient and proportional coefficient of auxiliary frequency controller, respectively; $T_{d}$ is the total time delay constant; and $H_{W T}$ is the mechanical inertia constant of the WT.

Equation (2) shows that various power signals of WT are coupled together through rotor speed and system frequency. Therefore, the output active power not only depends on the system frequency deviation, but also is related to the change of rotor speed. To analyze the influence of the change in rotor speed, power response curves of WT in the simulation with the change of rotor speed and time are shown in Fig. 2.

Fig. 2 Power response curves of WT. (a) P-ω_r. (b) $Δ P_{W e} - t$ .

In Fig. 2(a), when a frequency fault occurs, $P_{W m}$ and $P_{M P P T}$ decrease from the stable operating point A with the decrease of rotor speed. And the deviation of $P_{M P P T}$ $Δ P_{M P P T}$ counteracts the effect of $Δ P_{a d d}$ , which means that the change of the rotor speed operating point weakens the frequency support capability of WT. With the continuation of the process of rotor kinetic energy release, $Δ P_{M P P T}$ gradually plays a dominant role. Especially when the speed decreases to $ω_{r B}$ , $Δ P_{W e} = 0$ , and the WT cannot provide frequency support to the power grid. Until the rotor speed decreases to $ω_{r C}$ , the system reaches a new stable operating point C, where the deviations of $ω_{r}$ and $P_{W m}$ reach their maximum values.

In order to accurately evaluate the frequency support capability, it is necessary to reveal the relationship between the output active power of WF and the system frequency. Based on small-signal incremental analysis, the state space model of the rotor speed can be expressed as:

\{\begin{array}{l} s Δ X = A Δ X + B Δ U \\ Δ Y = C Δ X + D Δ U \\ A = \frac{\frac{1}{2} ρ S V^{3} \frac{\partial C_{P}}{\partial ω_{r_{}}} - 3 k_{o p t} ω_{r 0}^{2}}{2 H_{W T} ω_{r 0}} \\ B = \frac{K_{d} s + K_{p}}{2 H_{W T} ω_{r 0} (T_{d} s + 1)} \\ C = 3 k_{o p t} ω_{r 0}^{2} \\ D = - \frac{K_{d} s + K_{p}}{T_{d} s + 1} \end{array}

(3)

where $X = ω_{r}$ is the state variable; $U = f$ is the input variable; $Y = P_{W e}$ is the output variable; and $ω_{r 0}$ is the initial rotor speed of the WT before frequency fault.

As can be observed from Fig. 2(a), the difference between $P_{W e}$ and $P_{W m}$ gradually decreases, so $ω_{r}$ slows down and $Δ ω_{r}$ is small. In addition, when $Δ f$ reaches a steady state, $ω_{r}$ is usually still decreasing, so $ω_{r}$ does not decrease greatly in the time scale concerned. Accordingly, it can be observed that $ω_{r}$ is near the equilibrium, so the linearized model (3) is valid. Then, according to (3) and the transfer function $Δ Y = [C {(s - A)}^{- 1} B + D] Δ U$ , the relationship between $Δ P_{W e}$ and $Δ f$ during frequency regulation process can be obtained as:

\frac{Δ P_{W e} (s)}{Δ f (s)} = - \frac{2 H_{W T} ω_{r 0} s - \frac{1}{2} ρ S V^{3} \frac{\partial C_{P}}{\partial ω_{r_{}}}}{2 H_{W T} ω_{r 0} s - \frac{1}{2} ρ S V^{3} \frac{\partial C_{P}}{\partial ω_{r_{}}} + 3 k_{o p t} ω_{r 0}^{2}} \frac{K_{d} s + K_{p}}{T_{d} s + 1}

(4)

It can be concluded from (4) that the transfer function between $Δ P_{W e}$ and $Δ f$ is a complex second-order function considering the change in the operating point of WT instead of a simple differential and proportional function. The curve of $Δ P_{W e}$ obtained from the linearized model (4) is in good accordance to the curve considering the change in the operating point, i.e., the actual $Δ P_{W e}$ . Therefore, it is necessary to consider the change in operating state when analyzing the frequency response of the system including WF participating in frequency regulation, which can provide a basis for establishing an accurate frequency response model.

Combining (4) with ASF model, a full-order frequency response model for multi-machine power system considering the change in the operating point of WT can be obtained in Fig. 3, where $K_{m k}^{T}$ and $K_{m k}^{H}$ are the mechanical power gain coefficients of the $k^{t h}$ SG in the thermal power plants and hydropower plants, respectively; $R_{k}^{T}$ and $R_{k}^{H}$ are the droop coefficients of the speed governor of the $k^{t h}$ SG in the thermal power plants and hydropower plants, respectively; $T_{G k}^{T}$ and $T_{G k}^{H}$ are the time constants of the $k^{t h}$ SG prime mover and governor in the thermal power plants and hydropower plants, respectively; $F_{H k}^{T}$ is the ratio of power generated by the high-pressure turbine of the k^th SG in the thermal power plants; and $T_{R k}^{T}$ and $T_{W k}^{H}$ are the re-heater and water hammer effect time constants of the $k^{t h}$ SG in the thermal power plants and hydropower plants, respectively. Wind power, thermal power, and hydropower in Fig. 3 are examples, but the established full-order model is also applicable for other types of energy. The premise is that the relationship between the output active power and the frequency deviation can be expressed by transfer function. Therefore, compared with the classical single-machine equivalent model, the established full-order model has higher accuracy and better applicability.

Fig. 3 Full-order frequency response model.

The full-order model of the system frequency deviation is:

\begin{array}{l} Δ f (s) = - \frac{Δ P_{d}}{s} [2 H s + D + \frac{K_{m k}^{T} (1 + F_{H k}^{T} T_{R k}^{T} s)}{R_{k}^{T} (1 + T_{G k}^{T} s) (1 + T_{R k}^{T} s)} + \dots + \\ \frac{K_{m k}^{H} (1 - T_{W k}^{H} s)}{R_{k}^{H} (1 + T_{G k}^{H} s) (1 + 0.5 T_{W k}^{H} s)} + \dots + \\ \frac{2 H_{W T} ω_{r 0} s - \frac{1}{2} ρ S V^{3} \frac{\partial C_{P}}{\partial ω_{r}}}{2 H_{W T} ω_{r 0} s - \frac{1}{2} ρ S V^{3} \frac{\partial C_{P}}{\partial ω_{r}} + 3 k_{o p t} ω_{r 0}^{2}} \frac{K_{d} s + K_{p}}{T_{d} s + 1}] \end{array}

(5)

Considering the strong applicability of the established frequency response model, the full-order model of (5) can also be rewritten as a n-order model, where n depends on the number and type of generators connected to the grid:

Δ f (s) = - \frac{Δ P_{d}}{s} \frac{b_{n - 1} s^{n - 1} + \dots + b_{1} s + b_{0}}{a_{n} s^{n} + \dots + a_{2} s^{2} + a_{1} s + a_{0}}

(6)

where $a_{0}, a_{1}, \dots, a_{n}$ and $b_{0}, b_{1}, \dots, b_{n - 1}$ are the polynomial coefficients of the full-order model.

III. P-ROFR Model

A. Formulation of P-ROFR Model

The time-domain expression of the full-order frequency response model cannot be directly solved by inverse Laplace transform, thus it is challenging to explicitly express the frequency response indexes. The classical SFR model shows that the frequency response model during frequency faults can be expressed as a second-order model with two real poles or a pair of conjugate complex poles. Therefore, an ROFR model with second-order structure is proposed as:

Δ f_{R O} (s) = - \frac{Δ P_{d}}{s} \frac{d_{1} s + d_{0}}{s^{2} + c_{1} s + c_{0}}

(7)

where $c_{0}$ , $c_{1}$ , $d_{0}$ , and $d_{1}$ are the polynomial coefficients of the ROFR model.

Considering that different frequency response stages have different requirements for the ROFR model, the frequency response is divided into three stages, as shown in Fig. 4: the initial stage of frequency drop (called transient stage), the stage near FN (called intermediate stage), and the stage near $Δ f_{s s}$ (called steady-state stage). In order to obtain the frequency response model with high accuracy at each frequency response stage, different ROFR models are adopted in different stages. Therefore, a P-ROFR model $Δ f_{R O}^{P}$ composed of three ROFR models is proposed, as shown in (8). The ROFR models adopted in the transient stage $[0, t_{T}]$ , intermediate stage $[t_{T}, t_{S}]$ , and steady-state stage $[t_{S}, \infty)$ are called transient ROFR model $Δ f_{R O, T}$ , intermediate ROFR model $Δ f_{R O, I}$ , and steady-state ROFR model $Δ f_{R O, S}$ , respectively. The three ROFR models have the same second-order structure, but the coefficients $c_{0}, c_{1}, d_{0}, d_{1}$ of these models are different.

Fig. 4 P-ROFR model.

Δ f_{R O}^{P} (t) = \{\begin{matrix} \begin{matrix} L^{- 1} (Δ f_{R O, T}) = L^{- 1} (- \frac{Δ P_{d}}{s} \frac{d_{1, T} s + d_{0, T}}{s^{2} + c_{1, T} s + c_{0, T}}) & t \in [0, t_{T}] \end{matrix} \\ \begin{matrix} L^{- 1} (Δ f_{R O, I}) = L^{- 1} (- \frac{Δ P_{d}}{s} \frac{d_{1, I} s + d_{0, I}}{s^{2} + c_{1, I} s + c_{0, I}}) & t \in [t_{T}, t_{S}] \end{matrix} \\ \begin{matrix} L^{- 1} (Δ f_{R O, S}) = L^{- 1} (- \frac{Δ P_{d}}{s} \frac{d_{1, S} s + d_{0, S}}{s^{2} + c_{1, S} s + c_{0, S}}) & t \in [t_{S}, \infty) \end{matrix} \end{matrix}

(8)

where the piecewise point $t_{T}$ is the time when $L^{- 1} (Δ f_{R O, T} (s))$ intersects $L^{- 1} (Δ f_{R O, I} (s))$ ; and the piecewise point $t_{S}$ is the time when $L^{- 1} (Δ f_{R O, I} (s))$ intersects $L^{- 1} (Δ f_{R O, S} (s))$ .

In order to ensure that the P-ROFR model is as close as possible to the full-order frequency response model, a method to solve the coefficients of the P-ROFR model by minimizing the amplitude of the error model is proposed. Combining (6) and (7), the frequency-domain error model between the full-order frequency response model $Δ f$ and the ROFR model $Δ f_{R O}$ can be written as:

\begin{array}{l} δ (s) = Δ f (s) - Δ f_{R O} (s) = \\ - \frac{Δ P_{d}}{s} [\frac{(b_{n - 1} - a_{n} d_{1}) s^{n + 1} + (b_{n - 2} + b_{n - 1} c_{1} - a_{n - 1} d_{1} - a_{n} d_{0}) s^{n}}{(a_{n} s^{n} + \dots + a_{1} s + a_{0}) (s^{2} + c_{1} s + c_{0})} + \\ \frac{\sum_{i = 2}^{n - 1} (b_{i - 2} + b_{i - 1} c_{1} + b_{i} c_{0} - a_{i - 1} d_{1} - a_{i} d_{0}) s^{i}}{(a_{n} s^{n} + \dots + a_{1} s + a_{0}) (s^{2} + c_{1} s + c_{0})} + \\ \frac{(b_{0} c_{1} + b_{1} c_{0} - a_{0} d_{1} - a_{1} d_{0}) s + b_{0} c_{0} - a_{0} d_{0}}{(a_{n} s^{n} + \dots + a_{1} s + a_{0}) (s^{2} + c_{1} s + c_{0})}] \end{array}

(9)

It can be observed from (9) that if the coefficients of the numerator polynomial are close to 0, the frequency-domain characteristic of $Δ f_{R O}$ is similar to that of $Δ f$ . It means that the time-domain performance indexes of these two models should also be consistent. Therefore, in order to make the $n + 2$ coefficients of numerator polynomial in (9) be 0, $n + 2$ sets of equations are established, as shown in (10).

X θ = Y

(10)

where X, Y, and θ are the coefficient matrices, and θ can be obtained by the solution formula $θ = (X^{T} {X)}^{- 1} X^{T} Y$ of the least square method. X, θ, and Y can be expressed as:

X = [\begin{matrix} 0 & 0 & 0 & a_{n} \\ 0 & - b_{n - 1} & a_{n} & a_{n - 1} \\ - b_{n - 1} & - b_{n - 2} & a_{n - 1} & a_{n - 2} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ - b_{2} & - b_{1} & a_{2} & a_{1} \\ - b_{1} & - b_{0} & a_{1} & a_{0} \\ - b_{0} & 0 & a_{0} & 0 \end{matrix}]

(11)

θ = {[\begin{matrix} c_{0} & c_{1} & d_{0} & d_{1} \end{matrix}]}^{T}

(12)

Y = {[\begin{matrix} b_{n - 1} & b_{n - 2} & \dots & b_{0} & 0 & 0 \end{matrix}]}^{T}

(13)

According to the solution formula $θ = (X^{T} {X)}^{- 1} X^{T} Y$ , the coefficients of the intermediate ROFR model $θ_{I}$ are obtained by solving $n + 2$ sets of equations in the overdetermined equation (10), in order to minimize all the coefficients of the numerator in (9). The coefficients of the transient ROFR model $θ_{T}$ are obtained by solving the first four equations among the $n + 2$ sets of equations in (10), in order to make the coefficients of $s^{n + 2}$ , $s^{n + 1}$ , $s^{n}$ , and $s^{n - 1}$ of the numerator in (9) zero. The coefficients of the steady-state ROFR model $θ_{S}$ are obtained by solving the last four equations among the $n + 2$ sets of equations in (10), in order to make the coefficients of $s^{3}$ , $s^{2}$ , $s^{1}$ , and $s^{0}$ of the numerator in (9) zero.

The frequency-domain error model between the transient ROFR model and the full-order model can be expressed in (14). The frequency-domain error model between the steady-state ROFR model and the full-order model can be expressed in (15).

δ_{T} (s) = Δ f (s) - Δ f_{R O, T} (s) = - \frac{Δ P_{d}}{s} \frac{\sum_{i = 0}^{n - 3} M_{T, i} s^{i}}{\sum_{i = 0}^{n + 2} N_{T, i} s^{i}}

(14)

δ_{S} (s) = Δ f (s) - Δ f_{R O, S} (s) = - \frac{Δ P_{d}}{s} \frac{\sum_{i = 4}^{n + 1} M_{S, i} s^{i}}{\sum_{i = 0}^{n + 2} N_{S, i} s^{i}}

(15)

where $M_{T, i}$ , $N_{T, i}$ , $M_{S, i}$ , and $N_{S, i}$ are the coefficients of (10).

B. Mathematical Validation of P-ROFR Model

In this subsection, mathematical theorems are applied to prove that different frequency response stages should adopt the corresponding reduced-order models.

1)　Transient Stage

By performing inverse Laplace transform on (14), the time-domain error model between $Δ f_{R O, T}$ and $Δ f$ is described as:

δ_{T} (t) = \sum_{j = 1}^{n + 2} k_{j} e^{- μ_{j} t} + k_{0}

(16)

where $μ_{j}$ is the pole of (14); and $k_{j}$ $(j = 0,1, \dots, n + 2)$ is the amplitude of the signal.

In the transient stage, the frequency drops rapidly, so only the error model within a short time scale after the frequency drop occurs needs to be considered. The Taylor series expansion of (16) at $t = 0$ is given as:

δ_{T} (t) = \sum_{j = 0}^{n + 2} k_{j} + \sum_{j = 1}^{n + 2} (- k_{j} μ_{j} t + \frac{k_{j} μ_{j}^{2}}{2!} t^{2} - \frac{k_{j} μ_{j}^{3}}{3!} t^{3} + \frac{k_{j} μ_{j}^{4}}{4!} t^{4}) + o (t^{4})

(17)

where $o (t^{4})$ is the higher-order infinitesimal of $t^{4}$ .

Combining the initial value theorem of Laplace transform and (14), the initial value and the first four-order derivatives of $δ_{T} (t)$ are all 0 at $t = 0$ .

\{\begin{array}{l} δ_{T} (0) = \sum_{j = 0}^{n + 2} k_{j} = \underset{s \to \infty}{l i m} s δ_{T} (s) = 0 \\ \frac{d^{χ} δ_{T} (0)}{d t^{χ}} = \sum_{j = 1}^{n + 2} {(- 1)}^{χ} k_{j} μ_{j}^{χ} = \underset{s \to \infty}{l i m} s (s^{χ} δ_{T} (s)) = 0 \end{array}

(18)

where $χ$ is the order of the derivative ( $χ = 1, 2, 3, 4$ ).

Substituting (18) into (17), it can be observed that the error between $Δ f_{R O, T}$ and $Δ f$ is very small during a very short period after $t = 0_{+}$ .

δ_{T} (t) = o (t^{4})

(19)

Similarly, the error between $Δ f_{R O, I}$ and $Δ f$ , and the error between $Δ f_{R O, S}$ and $Δ f$ can be obtained respectively as:

δ_{I} (t) = o (t^{κ})

(20)

δ_{S} (t) = o (1)

(21)

where $o (t^{κ})$ is the higher-order infinitesimal of $t^{κ}$ ( $κ = 1, 2,$ $3, 4$ ); and $o (1)$ is the higher-order infinitesimal of 1.

Based on the higher-order infinitesimal theorem, it can be observed that $δ_{T} < < δ_{I}$ and $δ_{T} < < δ_{S}$ , so $Δ f_{R O, T}$ is closer to $Δ f$ than $Δ f_{R O, I}$ and $Δ f_{R O, S}$ in the transient stage.

2)　Intermediate Stage

It can be concluded in [

10] that the FN is related to the inertia of the system, the governor of the SG, the scheme of AFC in the WF, and the operating state of the unit, etc. Therefore, the reduced-order model in the intermediate stage should consider the effects of active power output from all devices. The process of solving

Δ f_{R O, T}

and

Δ f_{R O, S}

only uses four sets of equations in (10). However, the coefficients of the intermediate ROFR model

Δ f_{R O, I}

need to satisfy

n + 2

sets of equations in (10). Therefore,

Δ f_{R O, I}

contains all the information of

Δ f

, which is more suitable than

Δ f_{R O, T}

and

Δ f_{R O, S}

to express the frequency response characteristics in the intermediate stage.

3)　Steady-state Stage

Substituting $s = j ω$ into (9), (14), and (15), the ratios of the magnitude of $δ_{S}$ to that of $δ_{T}$ and that of $δ_{I}$ can be obtained, respectively, as:

|\frac{δ_{S}}{δ_{T}}| = ω^{4} \sum_{l = 0}^{n - 3} g (ω^{l})

(22)

|\frac{δ_{S}}{δ_{I}}| = ω^{κ} \sum_{l = 0}^{n - 3} h (ω^{l})

(23)

where $ω$ is the oscillation frequency in rad/s; and $g (\cdot)$ and $h (\cdot)$ are the polynomial functions of $ω$ .

In the steady-state stage, the oscillation frequency $ω < < 1$ , so $|δ_{S} / δ_{T}|$ and $|δ_{S} / δ_{I}|$ are both much less than 1. Therefore, $Δ f_{R O, S}$ is closer to $Δ f$ than $Δ f_{R O, T}$ and $Δ f_{R O, I}$ .

C. Numerical Validation of P-ROFR Model

The approximation degree of $Δ f_{R O, T}$ , $f_{R O, I}$ , $Δ f_{R O, S}$ , and P-ROFR model $Δ f_{R O}^{P}$ to $Δ f$ is compared in Fig. 5 and Table I, where $Δ f_{m a x}^{e r r o r}$ , $R o C o F_{m a x}^{e r r o r}$ , and $Δ f_{s s}^{e r r o r}$ are the error rates of $Δ f_{m a x}^{}$ , $R o C o F_{m a x}^{}$ , and $Δ f_{s s}^{}$ , respectively; and R² is an index characterizing the overall fitting degree. And the closer the value of R² is to 1, the better the fitting degree of the reduced-order model to the full-order model is. It can be observed that in the transient stage, compared with $Δ f_{R O, I}$ and $Δ f_{R O, S}$ , $Δ f_{R O, T}$ can more accurately characterize the frequency change trend after the frequency drop. In the intermediate stage, $Δ f_{R O, I}$ can represent the maximum frequency deviation ( $Δ f_{m a x}$ ) more accurately than $Δ f_{R O, T}$ and $Δ f_{R O, S}$ . In the steady-state stage, $Δ f_{s s}$ can be obtained by $Δ f_{R O, S}$ with higher accuracy than $Δ f_{R O, T}$ and $Δ f_{R O, I}$ . Moreover, the P-ROFR model composed of these three reduced-order models has the highest accuracy in each frequency response stage, and its R² is the closest to 1, which proves the effectiveness of the proposed P-ROFR model.

Fig. 5 Comparison of ROFR models. (a) Whole stage. (b) Transient stage. (c) Intermediate stage.

TABLE I SFR Results of Reduced-order Model

Model	$Δ f_{m a x}^{e r r o r}$ (%)	$R o C o F_{m a x}^{e r r o r}$ (%)	$Δ f_{s s}^{e r r o r}$ (%)	R²
$Δ f_{R O, T}$	-24.62	5.560×10^-10	1.2289×10²	-9.6938
$Δ f_{R O, I}$	-1.31	-1.789×10¹	-1.0480×10¹	0.9665
$Δ f_{R O, S}$	-4.04	-3.723×10¹	-7.5900×10^-10	0.9693
$Δ f_{R O}^{P}$	-1.31	5.560×10^-10	-7.5900×10^-10	0.9977

IV. O-AFC Scheme with Inner- and Outer-layerOptimization

Previous analysis shows that P-ROFR model is a second-order model, so the analytical expression of each frequency response index can be derived directly as a function of $K_{d}$ and $K_{p}$ , which provides the optimization objective for the O-AFC scheme of WF.

It is worth mentioning that when the additional power $Δ P_{a d d}$ is increased by enlarging $K_{d}$ and $K_{p}$ , the frequency support capability of WF is not necessarily improved. There are three main reasons.

Firstly, with the increase of $K_{d}$ and $K_{p}$ , more rotor kinetic energy is released, thus the output power of the MPPT controller has a stronger ability to weaken the frequency support capability of WF. Secondly, excessive $K_{d}$ and $K_{p}$ lead to the rotor speed exceeding the minimum speed limit, which will result in the WT off-grid operation. Finally, the time delay of the power response of WF nonlinearly affects the effect of $K_{d}$ and $K_{p}$ on the system frequency, and appropriate time delay can effectively improve the frequency response characteristics, as shown in Supplementary Material A. Therefore, in order to maximize the frequency support capability of the WF on the premise of ensuring the full release of rotor kinetic energy and the full use of time delay, an O-AFC scheme is proposed in this section.

A. Optimization Objective of O-AFC Scheme

Considering that the output active power of the WF mainly affects the RoCoF after a time delay, the RoCoF_avg is regarded as one of the optimization objectives of the O-AFC scheme in this paper. At the same time, since $Δ f_{m a x}$ and $Δ f_{s s}$ are significant indexes to evaluate the frequency stability of the system, the objective function is given as:

\begin{matrix} m i n J (K_{d}, K_{p}) = η_{1} \end{matrix} |R o C o F_{a v g}| + η_{2} |Δ f_{m a x}| + η_{3} |Δ f_{s s}|

(24)

where $η_{1}, η_{2}, η_{3} \in (0,1)$ are the weighting coefficients, which can be configured according to different control requirements on the premise of satisfying $η_{1} + η_{2} + η_{3} = 1$ . The time scale of RoCoF_avg is $t_{n a d i r} / 3$ , where t_nadir is the time to reach the FN. The explicit expressions of RoCoF_avg, $Δ f_{m a x}$ , and $Δ f_{s s}$ as functions of $K_{d}$ and $K_{p}$ are obtained from the time-domain expressions of $Δ f_{R O, T}$ , $Δ f_{R O, I}$ , and $Δ f_{R O, S}$ , respectively. The structure of the proposed $Δ f_{R O, T}$ , $Δ f_{R O, I}$ , and $Δ f_{R O, S}$ is consistent with the classical second-order SFR model. Therefore, based on the solution formula of the frequency response indexes in [

23], the analytical expression of each index in (24) can be derived as:

R o C o F_{a v g} = \frac{Δ f_{R O, T} (t_{n a d i r} / 3)}{t_{n a d i r} / 3} = \frac{- 3 Δ P_{d} (λ_{2, T} - λ_{1, T}) [- \frac{r_{1, T}}{λ_{1, T}} {(\frac{r_{1, T}}{r_{2, T}})}^{\frac{λ_{1, T}}{3 (λ_{2, T} - λ_{1, T})}} + \frac{r_{2, T}}{λ_{2, T}} {(\frac{r_{1, T}}{r_{2, T}})}^{\frac{λ_{2, T}}{3 (λ_{2, T} - λ_{1, T})}} + \frac{d_{0, T}}{λ_{1, T} λ_{2, T}}]}{l n r_{1, T} - l n r_{2, T}}

(25)

\begin{array}{l} Δ f_{m a x} = Δ f_{R O, I} (t_{n a d i r}) = \\ - Δ P_{d} [- \frac{r_{1, I}}{λ_{1, I}} {(\frac{r_{1, I}}{r_{2, I}})}^{\frac{λ_{1, I}}{λ_{2, I} - λ_{1, I}}} + \frac{r_{2, I}}{λ_{2, I}} {(\frac{r_{1, I}}{r_{2, I}})}^{\frac{λ_{2, I}}{λ_{2, I} - λ_{1, I}}} + \frac{d_{0, I}}{λ_{1, I} λ_{2, I}}] \end{array}

(26)

Δ f_{s s} = Δ f_{R O, S} (t_{e n d}) = - Δ P_{d} (- \frac{r_{1, S}}{λ_{1, S}} e^{λ_{1, S} t_{e n d}} + \frac{r_{2, S}}{λ_{2, S}} e^{λ_{2, S} t_{e n d}} + \frac{d_{0, S}}{λ_{1, S} λ_{2, S}})

(27)

where $t_{e n d}$ is the end time of frequency regulation, which is set to be 20 s; and the subscripts T, I, S represent the transient, intermediate, and steady-state stages, respectively. The detailed expressions of $λ_{1}$ , $λ_{2}$ , $r_{1}$ , and $r_{2}$ are given in Supplementary Material B, which are composed of the coefficients c₀, c₁, d₀, and d₁ of the ROFR model. It is worth mentioning that the analytical expressions of c₀, c₁, d₀, and $d_{1}$ are explicitly expressed as polynomial functions of $K_{d}$ and $K_{p}$ by solving (10). The expressions of $c (K_{d}, K_{p})$ and $d (K_{d}, K_{p})$ are slightly complicated, but their accuracy is high because all parameters of the full-order model are preserved during the order reduction to improve the accuracy. And since the polynomial function is elementary, so the running time and the computational burden are not increased.

B. Constraint of O-AFC Scheme

1) Rotor speed constraint: according to the analysis of Fig. 2, $ω_{r}$ gradually decreases to a steady state after a frequency fault occurs, so the deviation of $ω_{r, s s}$ is the largest. Therefore, in order to avoid the large linearization errors of $ω_{r, s s}$ derived from (3), the relationship between $ω_{r, s s}$ and control coefficient is derived based on the power balance, without linearization.

When the system reaches a steady state, $P_{W m}$ is equal to $P_{W e}$ .

P_{W m, s s} = P_{W e, s s}

(28)

where $P_{W m, s s}$ and $P_{W e, s s}$ are the steady states of $P_{W m}$ and $P_{W e}$ , respectively.

When the pitch angle and wind speed are known, P_Wm can be approximated by a second-order polynomial of $ω_{r}$ [

9] as:

P_{W m} = k_{m 2} ω_{r}^{2} + k_{m 1} ω_{r} + k_{m 0}

(29)

where $k_{m 0}$ , $k_{m 1}$ , and $k_{m 2}$ are the constant coefficients.

When the system is in a steady state, the differential module of the auxiliary frequency controller cannot support the frequency, so $P_{W e, s s}$ can be expressed as:

P_{W e, s s} = k_{o p t} ω_{r, s s}^{3} - K_{p} Δ f_{s s}

(30)

Combining (28), (29), and (30), $Δ f_{s s}$ can be obtained as:

Δ f_{s s} = \frac{k_{o p t} ω_{r, s s}^{3} - k_{m 2} ω_{r, s s}^{2} - k_{m 1} ω_{r, s s} - k_{m 0}}{K_{p}}

(31)

According to the SFR model (5), $Δ f_{s s}$ can also be expressed as:

Δ f_{s s} = \frac{- Δ P_{d} + k_{o p t} ω_{r, s s}^{3} - k_{o p t} ω_{r 0}^{3}}{D + \frac{K_{m 1}}{R_{1}} + \dots + \frac{K_{m k}}{R_{k}} + K_{p}}

(32)

Combining (31) and (32), the relationship between $ω_{r, s s}$ and $K_{p}$ is expressed as:

K_{p} = \frac{(k_{o p t} ω_{r, s s}^{3} - k_{m 2} ω_{r, s s}^{2} - k_{m 1} ω_{r, s s} - k_{m 0}) (D + \frac{K_{m 1}}{R_{1}} + \dots + \frac{K_{m k}}{R_{k}})}{- Δ P_{d} + k_{m 2} ω_{r, s s}^{2} + k_{m 1} ω_{r, s s} + k_{m 0} - k_{o p t} ω_{r 0}^{3}}

(33)

According to (33), the relationship between $ω_{r, s s}$ and $K_{p}$ can be shown in Supplementary Material C. It can be observed that $ω_{r, s s}$ gradually decreases with the increase of $K_{p}$ . Therefore, the constraint that $ω_{r, s s}$ should be greater than $ω_{r m i n}$ can be converted to a constraint on $K_{p}$ , as shown in (34).

K_{p} \leq \frac{(k_{o p t} ω_{r m i n}^{3} - k_{m 2} ω_{r m i n}^{2} - k_{m 1} ω_{r m i n} - k_{m 0}) (D + \frac{K_{m 1}}{R_{1}} + \dots + \frac{K_{m k}}{R_{k}})}{- Δ P_{d} + k_{m 2} ω_{r m i n}^{2} + k_{m 1} ω_{r m i n} + k_{m 0} - k_{o p t} ω_{r 0}^{3}}

(34)

where $ω_{r m i n}$ is the lower limit of the rotor speed.

It is worth mentioning that the linearization is not included in the derivation of (28)-(34). The nonlinear characteristics of $ω_{r}$ are preserved in (34), so there is no linearization error in the constraint of $ω_{r, s s}$ .

2) Frequency response index constraint: considering that the weighting coefficients in (24) are determined by the control requirements, excessive attention of the system operator to a certain frequency response index and setting of a large weighting coefficient may cause other indexes to exceed the limit. Therefore, in order to ensure the frequency operating within a safe range, it is necessary to satisfy the maximum frequency deviation constraint and the steady-state frequency deviation constraint:

|Δ f_{m a x}| < Δ f_{m a x, l i m i t}

(35)

|Δ f_{s s}| < Δ f_{s s, l i m i t}

(36)

where $Δ f_{m a x, l i m i t}$ is the maximum frequency deviation limit;and $Δ f_{s s, l i m i t}$ is the steady-state frequency deviation limit.

C. Implementation of O-AFC Scheme

Considering that the effect of time delay is not as great as that of $K_{d}$ and $K_{p}$ , an O-AFC scheme with two-layer optimization is proposed, in which the inner layer optimizes the control coefficients $K_{d}$ and $K_{p}$ , and the outer layer optimizes the time delay $T_{d}$ .

The detailed implementation of the proposed inner- and outer-layer optimizations is shown in Supplementary Material D. When a frequency fault occurs, $T_{d}$ is increased by $Δ T$ each time in the outer optimization and the changed $T_{d}$ is substituted into the optimization as a known quantity. In the inner-layer optimization, particle swarm optimization algorithm is adopted to minimize the objective function (24) under the constraints (34)-(36), so as to obtain the optimized $K_{d}$ and $K_{p}$ under this $T_{d}$ . Finally, the optimized $T_{d}$ , $K_{d}$ , and $K_{p}$ can be obtained by comparing the minimum objective function under different $T_{d}$ .

V. Case Study

Two test power systems based on a two-area four-machine test system are built in MATLAB/Simulink, as shown in Fig. 6. In test system 1, G1-G4 are thermal generators with a total rated capacity of 2800 MW, and a WF is connected to bus B2, which has 469 doubly-fed induction generators with rated capacity of 1.5 MW. In test system 2, thermal generator G4 is replaced by a hydroelectric generator with the same rated capacity. The frequency base is 60 Hz and the power disturbance is set to be 270 MW.

Fig. 6 Configuration of two-area four-machine test system.

A. Verification of P-ROFR Model in Test System 1

In order to verify the accuracy and parameter applicability of the proposed P-ROFR model, two model parameter distribution cases of the SG are verified in test system 1, namely, Case A and Case B. The model parameters in Case A are uniformly distributed, as shown in Table II. The model parameters in Case B are set as polarized distribution, i.e., parameters are generally large or small values, and rarely intermediate values, as shown in Table III.

TABLE II Uniform Distribution Parameters in Case A

SG	$K_{m}^{T}$	$R_{}^{T}$	$F_{H}^{T}$	$T_{R}^{T}$	$T_{G}^{T}$
G1	0.15	0.10	0.39	9.10	0.21
G2	0.21	0.05	0.29	12.20	0.17
G3	0.29	0.08	0.25	6.30	0.27
G4	0.35	0.03	0.17	14.00	0.24

TABLE III Polarization Distribution Parameters in Case B

SG	$K_{m}^{T}$	$R_{}^{T}$	$F_{H}^{T}$	$T_{R}^{T}$	$T_{G}^{T}$
G1	0.16	0.10	0.39	7.80	0.30
G2	0.18	0.03	0.29	13.60	0.15
G3	0.32	0.04	0.15	14.20	0.28
G4	0.34	0.08	0.17	9.10	0.17

1)　Case A

The frequency response curves obtained from the full-order frequency response model, the improved SFR (I-SFR) model of [

15], and the proposed P-ROFR model in Case A are depicted in Fig. 7. Table IV shows the frequency response results under Case A. It can be observed that the

R o C o F_{a v g}^{e r r o r}

of the I-SFR model gradually increases during frequency fault in Fig. 7, and

Δ f_{m a x}^{e r r o r}_{}

between I-SFR model and full-order model exceeds 3% in Table IV, which are both caused by oversimplification of I-SFR modeling.

Fig. 7 SFR in Case A.

TABLE IV SFR Results in Case A

Model	$Δ f_{m a x}^{e r r o r}$ (%)	$R o C o F_{a v g}^{e r r o r}$ (%)	$Δ f_{s s}^{e r r o r}$ (%)	R²
I-SFR	-3.24	-1.51	-0.44	0.9751
P-ROFR	8.46×10^-3	-2.31×10^-3	-0.14	0.9974

However, the fitting degree between the proposed P-ROFR model and the full-order model is high throughout the frequency response process, and the error of each frequency response index is small enough to be ignored. Therefore, compared with I-SFR model, P-ROFR model not only has the same simplicity as the second-order model, but also exhibits the high accuracy in each frequency response stage.

2)　Case B

The SFR in Case B is shown in Fig. 8, and the SFR results in Case B are shown in Table V. It can be observed that the proposed P-ROFR model can better predict the frequency response characteristics of the power system. In particular, $Δ f_{m a x}^{e r r o r}$ between P-ROFR model and full-order model is much smaller than that between I-SFR model and full-order model, indicating that the P-ROFR model can accurately predict the FN. In addition, compared with Table IV, the accuracy of the P-ROFR model in Table V is still higher. Compared with Table IV, the accuracy of I-SFR model in Table V is greatly affected, especially the $R o C o F_{a v g}^{e r r o r}$ , which is mainly caused by oversimplification such as ignoring the large governor time constant of SG in Case B. The above analysis further shows that the accuracy of the I-SFR model is easily affected by the uneven distribution of parameters. In contrast, the proposed P-ROFR model has strong parameter applicability in different cases and has high accuracy during the whole frequency response process.

Fig. 8 SFR under Case B.

TABLE V SFR Results in Case B

Model	$Δ f_{m a x}^{e r r o r}$ (%)	$R o C o F_{a v g}^{e r r o r}$ (%)	$Δ f_{s s}^{e r r o r}$ (%)	R²
I-SFR	-7.0300	-4.4000	-0.45	0.9396
P-ROFR	0.0763	-0.0171	-0.14	0.9965

To further verify the prediction performance of the proposed P-ROFR model for FN, the relationships between $Δ f_{m a x}^{e r r o r}$ and K_d, K_p in Case A and Case B are shown in Figs. 9 and 10, respectively. It can be observed that in Case A and Case B, $Δ f_{m a x}^{e r r o r}$ between the P-ROFR model and the full-order model is significantly smaller than that between the I-SFR model and the full-order model. It indicates that the proposed P-ROFR model can achieve the high precision prediction of FN.

Fig. 9 Relationship between $Δ f_{m a x}^{e r r o r}$ and K_d, K_p under Case A. (a) P-ROFR model versus full-order model. (b) I-SFR model versus full-order model.

Fig. 10 Relationship between $Δ f_{m a x}^{e r r o r}$ and K_d, K_p in Case B. (a) P-ROFR model versus full-order model. (b) I-SFR model versus full-order model.

B. Verification of P-ROFR Model in Test System 2

To verify the simplification, accuracy, model applicability, and practical feasibility of the proposed P-ROFR model, the dominant pole retention method (DPRM) [

24] and the proposed method are used to reduce the order of the full-order model.

Frequency response models before and after order reduction are listed in Table VI. It can be observed that there may be fewer non-dominant poles in the full-order frequency response model, so it is very difficult to obtain a second-order model by using the DPRM. In contrast, when the piecewise reduced-order method is applied to the full-order model, the frequency response model realizes the order reduction from the 11^th order to the 2^nd order, which verifies the high simplification ability of the proposed P-ROFR model.

TABLE VI Frequency Response Models Before and After Order Reduction

Model

Expression

Full-order frequency

response model

Δ f (t) = L^{- 1} (\frac{- 0.00845 \prod_{i = 1}^{10} (s + z_{i})}{s \prod_{i = 1}^{11} (s + p_{i})})

DPRM reduced-order model

Δ f (t) = L^{- 1} (\frac{- 0.00845 \prod_{i = 1}^{7} (s + z_{i})}{s \prod_{i = 1}^{8} (s + p_{i})})

P-ROFR model

Δ f (t) = \{\begin{array}{l} L^{- 1} (- \frac{0.5 s + 10.40}{s (s^{2} + 21.55 s + 47.46)}) 0 < t \leq 0.445 \\ L^{- 1} (- \frac{0.49 s + 0.064}{s (s^{2} + 2.08 s + 0.45)}) 0.445 < t \leq 6.282 \\ L^{- 1} (- \frac{0.44 s + 0.061}{s (s^{2} + 1.83 s + 0.42)}) 6.282 < t < \infty \end{array}

Note: z_i and p_i are zeros and poles of frequency-domain expression, respectively.

The comparative results of simulation model, DPRM reduced-order model, and the proposed P-ROFR model are shown in Fig. 11. It can be observed from Fig. 11 that ignoring non-dominant poles greatly weakens the accuracy of the DPRM reduced-order model and causes large steady-state error. However, the proposed P-ROFR model has little difference with the actual simulation model in test system 2 with a variety of heterogeneous devices, which indicates that the proposed P-ROFR model has satisfactory accuracy and high model applicability. In addition, comparing Fig. 11(a) with Fig. 11(b), it can be observed that the accuracy of the proposed P-ROFR model is still higher under different operating conditions, which further proves the practical feasibility of the proposed method.

Fig. 11 SFR under different operation conditions. (a) $V = 9$ m/s and $P_{d} = 270$ MW. (b) $V = 12$ m/s and $P_{d} = 380$ MW.

Moreover, based on the simplicity of the second-order model, the analytical expression of the frequency indexes can be derived from the proposed P-ROFR model, which provides a basis for the O-AFC scheme of WF.

C. Verification of O-AFC Scheme in Test System 2

The coefficients of O-AFC scheme are obtained by two-layer optimization on the premise of ensuring the full release of the rotor kinetic energy and the full use of the effect of time delay. In the inner-layer optimization, the weighting coefficients of the optimization objective are $η_{1} = 0.3$ , $η_{2} = 0.6$ , and $η_{3} = 0.1$ , $ω_{r m i n}$ is set to be 0.85 p.u., and $Δ f_{m a x, l i m i t}$ and $Δ f_{s s, l i m i t}$ are set to be 0.5 Hz and 0.2 Hz, respectively. In the outer-layer optimization, $Δ T$ is set to be 0.05 s. In order to verify the effectiveness of O-AFC scheme, scheme 1, scheme 2, and O-AFC scheme are applied in test system 2 in two test cases, namely Case C and Case D. Scheme 1 is based on the improved frequency regulation scheme of [

23]. Scheme 2 ignores the effect of time delay on frequency regulation based on O-AFC scheme, which means that

K_{d}

and

K_{p}

of scheme 2 are the same as that of O-AFC scheme, but

T_{d}

of scheme 2 is inherent delay

T_{d 0}

(0.10 s). The optimized coefficients of three schemes in Case C and Case D (as indicated by subscripts C and D, respectively) are shown in Table VII.

TABLE VII Frequency Regulation Scheme in Case C and Case D

Scheme	K_d,C	K_p,C	T_d,C (s)	K_d,D	K_p,D	T_d,D (s)
Scheme 1	17.0	9.1	0.10	31.2	18.7	0.10
Scheme 2	37.1	15.8	0.10	38.5	45.2	0.10
O-AFC scheme	37.1	15.8	1.35	38.5	45.2	0.90

Note: in Case C, $V = 9$ m/s, $Δ P_{d} = 270$ MW, and in Case D, $V = 10$ m/s, $Δ P_{d} = 220$ MW.

1)　Case C

$K_{d, C}$ , $K_{p, C}$ , and $T_{d, C}$ in Case C in Table VII are applied to the test system 2. The curves of SFR, rotor speed, and output active power of WF are shown in Fig. 12.

Fig. 12 System performance in Case C ( $V = 9$ m/s, $Δ P_{d} = 270$ MW). (a) SFR. (b) Rotor speed of aggregated WF model. (c) Power response of WF.

It can be observed that the rotor kinetic energy of scheme 1 is not fully released, and the output active power is not enough to maximize the frequency support capability of the WF. This is because scheme 1 ignores the negative effect of the change in the operating point of WT on $Δ P_{a d d}$ , and the designed $K_{d}$ and $K_{p}$ are based on $Δ P_{a d d}$ instead of $Δ P_{W e}$ . Compared with scheme 1, the proposed O-AFC scheme has better frequency response characteristics by fully releasing rotor kinetic energy, especially the FN is significantly improved.

It can be observed from Fig. 12 that even though the optimized $K_{d}$ and $K_{p}$ of scheme 2 are the same as those of the O-AFC scheme, the frequency support capability of the O-AFC scheme with a larger time delay is better. And it can be observed from Fig. 12(c) that the maximum active power output by WF in O-AFC scheme is 0.192 p.u., which is 16.35% greater than that in scheme 2. This is because the time delay causes the frequency to drop faster in the initial time, so $Δ P_{a d d}$ increases according to the equation $Δ P_{a d d} = - K_{d} d Δ f / d t - K_{p} Δ f$ . And it can also be observed from Fig. 12(a) that the RoCoF of the O-AFC scheme suddenly slows down at 1.35 s, and the FN is increased compared with scheme 2. This is because the time delay is equivalent to delaying the output active power of WF, so the frequency support capability of the O-AFC scheme near the FN is greater than that of scheme 2. Therefore, the O-AFC scheme considering the effect of time delay has better frequency regulation performance.

Moreover, the optimization objective of the O-AFC scheme is derived from the proposed P-ROFR model, so the effectiveness of the O-AFC scheme also proves the high accuracy and practical feasibility of the P-ROFR model.

2)　Case D

$K_{d, D}$ , $K_{p, D}$ , and $T_{d, D}$ in Case D in Table VII are also applied to test system 2. The system performance in Case D is shown in Fig. 13. It can be observed that the frequency performance of the proposed O-AFC scheme is still better than that of scheme 1 and scheme 2 under different operating conditions and power disturbances.

Fig. 13 System performance in Case D ( $V = 10$ m/s, $Δ P_{d} = 220$ MW). (a) SFR. (b) Rotor speed of aggregated WF model. (c) Power response of WF.

And it can be observed from Table VII, Fig. 12, and Fig. 13 that when the wind speed is larger and the power disturbance is smaller, the optimized control coefficients $K_{d}$ and $K_{p}$ are larger. This is because the initial rotor speed of WT increases with the increase of wind speed, so when the rotor speed drops to the lower limit ( $ω_{r m i n} = 0.85$ p.u.), the deviation of the rotor speed is relatively large, as shown in Fig. 12(b) and Fig. 13(b). This means that in order to fully release the rotor kinetic energy, larger $K_{d}$ and $K_{p}$ are required, which effectively maximize the frequency support capability of the WF. In addition, it can be observed from (5) that the system frequency drops less when the power disturbance is smaller, so the output active power by releasing the rotor kinetic energy is reduced. Therefore, in order to fully release the rotor kinetic energy, the optimized $K_{d}$ and $K_{p}$ should be larger under the smaller power disturbance. To sum up, the optimized $K_{d}$ and $K_{p}$ will change with the operating conditions of WT and the power disturbances of the system, in order to ensure the full release of rotor kinetic energy.

To sum up, compared with I-SFR model, the proposed P-ROFR model achieves higher accuracy and higher parameter adaptability, especially higher-precision prediction of the FN. Compared with the DPRM reduced-order model, the proposed P-ROFR model not only achieves satisfactory accuracy but also has a simpler second-order structure, which can derive the explicit expression of frequency response indexes for the O-AFC scheme. Compared with scheme 1 and scheme 2, the proposed O-AFC scheme based on P-ROFR model effectively improves the frequency support capability of the WF while ensuring the full release of the rotor kinetic energy and the safe operation of the system. And the proposed O-AFC scheme makes full use of the effect of time delay on increasing the active power output by the WF and delaying the support time of the FN.

VI. Conclusion

In this paper, a P-ROFR model and an O-AFC scheme based on the proposed P-ROFR model have been proposed. The negative effect of the change in the operating point of WT on the frequency support capability has been analyzed, and a full-order frequency response model with high applicability has been established. A ROFR model with second-order structure has been proposed, thus the order of full-order model has been reduced. A P-ROFR model composed of three ROFR models has been further proposed, thus the accuracy has been improved throughout the frequency response process. Based on the P-ROFR model, the relationship between each frequency response index and the coefficients of auxiliary frequency controller has been explicitly expressed, which has provided the optimization objective for the O-AFC scheme. Then, the O-AFC scheme with inner- and outer-layer optimizations has been proposed, thus the frequency support capability of WF has been improved on the premise of ensuring the full release of the rotor kinetic energy and the full use of the effect of time delay on frequency regulation. Finally, the high accuracy, simplicity, applicability, and practical feasibility of the P-ROFR model and the effectiveness of the O-AFC scheme have been verified by the simulation studies.

References

P. Bao, W. Zhang, and Y. Zhang, “Secondary frequency control considering optimized power support from virtual power plant containing aluminum smelter loads through VSC-HVDC link,” Journal of Modern Power Systems and Clean Energy, vol. 11, no. 1, pp. 355-367, Jan. 2023. [Baidu Scholar]

J. Khazaei, Z. Tu, and W. Liu, “Small-signal modeling and analysis of virtual inertia-based PV systems,” IEEE Transactions on Energy Conversion, vol. 35, no. 2, pp. 1129-1138, Jun. 2020. [Baidu Scholar]

H. Zhang, W. Xiang, W. Lin et al., “Grid forming converters in renewable energy sources dominated power grid: control strategy, stability, application, and challenges,” Journal of Modern Power Systems and Clean Energy, vol. 9, no. 6, pp. 1239-1256, Nov. 2021. [Baidu Scholar]

S. Ghosh, S. Kamalasadan, N. Senroy et al., “Doubly-fed induction generator (DFIG)-based wind farm control framework for primary frequency and inertial response application,” IEEE Transactions on Power Systems, vol. 31, no. 3, pp. 1861-1871, May 2016. [Baidu Scholar]

C. Zhao, D. Sun, L. Zhang et al., “Two-stage optimal re/active power control of large-scale wind farm to mitigate voltage dip induced frequency excursion,” IEEE Transactions on Sustainable Energy, vol. 14, no. 1, pp. 730-733, Jan. 2023. [Baidu Scholar]

M. Toulabi, S. Bahrami, and A. M. Ranjbar, “An input-to-state stability approach to inertial frequency response analysis of doubly-fed induction generator-based wind turbines,” IEEE Transactions on Energy Conversion, vol. 32, no. 4, pp. 1418-1431, Dec. 2017. [Baidu Scholar]

S. Wang and K. Tomsovic, “Fast frequency support from wind turbine generators with auxiliary dynamic demand control,” IEEE Transactions on Power Systems, vol. 34, no. 5, pp. 3340-3348, Sept. 2019. [Baidu Scholar]

R. Yang, G. Shi, X. Cai et al., “Autonomous synchronizing and frequency response control of multi-terminal DC systems with wind farm integration,” IEEE Transactions on Sustainable Energy, vol. 11, no. 4, pp. 2504-2514, Oct. 2020. [Baidu Scholar]

H. Ye, W. Pei, and Z. Qi, “Analytical modeling of inertial and droop responses from a wind farm for short-term frequency regulation in power systems,” IEEE Transactions on Power Systems, vol. 31, no. 5, pp. 3414-3423, Sept. 2016. [Baidu Scholar]

L. Liu, W. Li, Y. Ba et al., “An analytical model for frequency nadir prediction following a major disturbance,” IEEE Transactions on Power Systems, vol. 35, no. 4, pp. 2527-2536, Jul. 2020. [Baidu Scholar]

M. Sun, Y. Min, X. Xiong et al., “Practical realization of optimal auxiliary frequency control strategy of wind turbine generator,” Journal of Modern Power Systems and Clean Energy, vol. 10, no. 3, pp. 617-626, May 2022. [Baidu Scholar]

M. L. Chan, R. D. Dunlop, and F. Schweppe, “Dynamic equivalents for average system frequency behavior following major distribances,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-91, no. 4, pp. 1637-1642, Jul. 1972. [Baidu Scholar]

P. M. Anderson and M. Mirheydar, “A low-order system frequency response model,” IEEE Transactions on Power Systems, vol. 5, no. 3, pp. 720-729, Aug. 1990. [Baidu Scholar]

M. Sun, Y. Min, L. Chen et al., “Optimal auxiliary frequency control of wind turbine generators and coordination with synchronous generators,” CSEE Journal of Power and Energy Systems, vol. 7, no. 1, pp. 78-85, Jan. 2021. [Baidu Scholar]

M. Qing, F. Tang, F. Liu et al., “An analytical method for estimating the maximum penetration of DFIG considering frequency stability,” Sustainability, vol. 12, no. 23, pp. 1-21, Dec. 2020. [Baidu Scholar]

Y. Chen, W. Xu, Y. Liu et al., “Reduced-order system frequency response modeling for the power grid integrated with the type-II doubly-fed variable speed pumped storage units,” IEEE Transactions on Power Electronics, vol. 37, no. 9, pp. 10994-11006, Sept. 2022. [Baidu Scholar]

H. Huang, P. Ju, Y. Jin et al., “Generic system frequency response model for power grids with different generations,” IEEE Access, vol. 8, pp. 14314-14321, Feb. 2020. [Baidu Scholar]

C. Zhang, L. Liu, and H. Cheng et al., “Frequency-constrained co-planning of generation and energy storage with high-penetration renewable energy,” Journal of Modern Power Systems and Clean Energy, vol. 9, no. 4, pp. 760-775, Jul. 2021. [Baidu Scholar]

J. Huang, Z. Yang, J. Yu et al., “Frequency dynamics-constrained parameter design for fast frequency controller of wind turbine,” IEEE Transactions on Sustainable Energy, vol. 13, no. 1, pp. 31-43, Jan. 2022. [Baidu Scholar]

X. Ge, X. Zhu, Y. Fu et al., “Optimization of reserve with different time scales for wind-thermal power optimal scheduling considering dynamic deloading of wind turbines,” IEEE Transactions on Sustainable Energy, vol. 13, no. 4, pp. 2041-2050, Oct. 2022. [Baidu Scholar]

J. Huang, Z. Yang, J. Yu et al., “Optimization for DFIG fast frequency response with small-signal stability constraint,” IEEE Transactions on Energy Conversion, vol. 36, no. 3, pp. 2452-2462, Sept. 2021. [Baidu Scholar]

M. Sun, Y. Sun, L. Chen et al., “Novel temporary frequency support control strategy of wind turbine generator considering coordination with synchronous generator,” IEEE Transactions on Sustainable Energy, vol. 13, no. 2, pp. 1011-1020, Apr. 2022. [Baidu Scholar]

C. Zhao, D. Sun, X. Zhang et al., “A two-stage power distribution scheme of multiple wind farms participating in primary frequency regulation,” IEEE Transactions on Power Systems, doi: 10.1109/TPWRS.2022.3224971 [Baidu Scholar]

S. K. Tiwari and G. Kaur, “Improved reduced-order modeling using clustering method with dominant pole retention,” IETE Journal of Research, vol. 66, no. 1, pp. 42-52, Jan. 2020. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher