According to aerodynamic principles (Shkara et al., 2018), the power captured by WT can be expressed as P in = 1 2 ρ S v 3 = 1 2 ρ π R 2 v 3 (1) where ρ , S , and v represent air density, the area swept by the blade of WT, and wind speed, respectively; R is the radius of the blade of WT.

Furthermore, the output power of WT can be calculated as P out = P in ⋅ C p ( λ ,   θ ) = 1 2 ρ π R 2 v 3 C p ( λ ,   θ ) (2) where λ and θ denote tip speed ratio and the pitch angle of WT, respectively; λ = ω R / v , ω stands for the rotational angular speed of WT in rad/s; C p ( λ ,   θ ) is wind energy utilization coefficient (Slootweg et al., 2003), which satisfies the following expression: C p ( λ ,   θ ) = c 1 ⋅ ( c 2 λ i − θ c 3 − c 4 ) ⋅ e c 5 λ i + c 6 (3) 1 λ i = 1 λ + 0.08 β − 0.035 θ 3 + 1 (4) where c 1 = 0.5176 ; c 2 = 116 ; c 3 = 0.4 ; c 4 = 5 ; c 5 = − 21 ; c 6 = 0.0068 .

From Eqs. 3, 4, one can know that the wind energy utilization coefficient C p of WT is determined by λ and θ , and the relationships of C p between λ and θ are illustrated in Figure 1. It can be seen from Figure 1 that the wind energy utilization coefficient increases and then decreases with the increase of tip speed ratio for a certain pitch angle. Besides, the maximum wind energy utilization coefficient decreases with the increase of pitch angle. Therefore, maximum power output can be achieved by controlling the angular speed of WT under different wind speeds. Furthermore, the pitch angle θ of WT remains as 0° for capturing maximum power output, and rated wind speed is settled as 12 m/s in this paper.

As shown in Figure 2, a three-area four-terminal voltage source converter-based multi-terminal high voltage direct current system (VSC-MT-HVDCs) combined with WFs and AC system is adopted as the test system, which consists of two WFs, a VSC-MT-HVDCS, and an AC system (Vennelaganti and Chaudhuri, 2018; Xiong et al., 2021), among which each WF includes five equivalent WTs that transmit electricity to the AC system together. In this paper, the rotor angular speeds of WT1–WT5 in WF1 are determined as 0.90 pu, 0.95 pu, 0.85 pu, 1.00 pu, and 1.05 pu, respectively. Besides, the rotor angular speeds of WT6–WT10 are also settled as 0.90 pu, 0.95 pu, 0.85 pu, 1.00 pu, and 1.05 pu, respectively. Note that only WTs of WF1 participate in frequency support through step start-up adaptive inertial droop control scheme after frequency events (e.g., load increase or decrease) while all WTs of WF2 merely operate in maximum power point (MPP) states. Meanwhile, WTs in the WF1 are categorized into different clusters in a certain ratio according to their rotor angular speeds. For instance, under load increase, WT2, WT4, and WT5 belong to cluster 1 while cluster 2 involves WT1 and WT3 if the cluster classification ratio is settled as 6:4. On the contrary, under load decrease, WT1–WT3 are classified as cluster 1 while WT4 and WT5 are cluster 2 if the classification ratio is also 6:4.

For the ith WT, furthermore, the step start-up adaptive inertial droop control scheme is demonstrated in Figure 3, among which the active power reference value P ref , i of the ith WT mainly consists of three parts, i.e., MPP under real-time rotor angular speed ω i of this WT, inertial control, and droop control, which is reacted after a certain time delay Δ t for matching with WTs of another cluster. In other words, cluster 1 and cluster 2 are utilized to respectively support system frequency in a certain order. Meanwhile, both coefficients of inertial control and droop control are adaptively varied with ω i , respectively. Besides, f AC stands for the frequency of the AC system, which is measured in VSC1 in this paper; f ref = 50 Hz , ω min = 0.70 pu , and ω max = 1.20 pu .

Compared with the conventional droop control method that adopts constant coefficients, adaptive inertial droop control coefficients vary with real-time rotor angular speed ω i of WT, which can not only guarantee rotor operating within a safe condition but also effectively utilize wind energy and timely adjust states (i.e., absorb or generate power) of WT, among which the variation rate of AC system frequency f AC is used as the input signal of inertial control while the input signal of droop control is the frequency deviation. Owing to the frequency deviation being small while frequency variation rate being large in the primary stage of frequency response, inertial control can rapidly provide frequency support (Kayikci and Milanovic, 2009). Besides, droop control is utilized to simulate the primary frequency regulation capability of SGs, which can improve steady-state frequency in the later stage after frequency events.

Additionally, more information about modeling parameters of the whole system and control strategy of VSC stations can be obtained from the literature (Xiong et al., 2021).

Inspired by the courtship behaviors of spiders, Hayyolalam and PourhajiKazem (2020) proposed a novel meta-algorithm, namely, BWOA, which obtained satisfactory performance not only in test functions but also in engineering optimization applications (Adrián et al., 2020).

In order to ensure the global searching ability of BWOA, the population with N individuals is initialized as X = [ x 1,1 x 2,1 ⋮ x l , 1 ⋮ x N , 1 x 1,2 x 2,2 ⋱ x l , 2 ⋱ x N , 2 ⋯ x 1 , D ⋯ x 2 , D ⋱ ⋮ ⋯ x l , D ⋱ ⋮ ⋯ x N , D ] (1) where D represents the dimension of the problem to be solved; each row of candidate solution matrix X determines a current location (corresponding to a solution) of a spider. Thus, the lth individual can be presented as X l → = ( U B − L B ) ⋅ r a n d + L B = [ x l , 1 , x l , 2 , ⋯ x l , D ] (2) where U B and L B   denote the upper bound vector and the lower bound vector of optimized parameters; r a n d is equal to a random number from 0 to 1.

Besides, the fundamental update rules of BWOA during iteration can be mathematically described as { X l → ( t + 1 ) = X b → ( t ) − m ⋅ X r → ( t )         i f   r a n d ( ) ≤ 0.3 X l → ( t + 1 ) = X b → ( t ) − c o s ( 2 π β ) ⋅ X l → ( t )       o t h e r s (3) where X l → ( t + 1 ) and X l → ( t ) denote new solution and the old one for the lth individual in the ( t + 1 ) th iteration, respectively; X b → ( t ) and X l → ( t ) stand for the best solution and a random solution in the previous iteration, respectively, with r ≠ l ; m is defined as a random float number from 0.4 to 0.9 while β is randomly generated in the interval of [−1.0,1.0].

Furthermore, pheromones of female spiders decide their mating rates with males, which can be calculated by p h ( l ) = F i t n e s s max − F i t n e s s ( l ) F i t n e s s max − F i t n e s s min (4) where p h ( l ) represents the pheromone value of the lth female spider, which is a float number from 0 to 1; F i t n e s s max , F i t n e s s min , and F i t n e s s ( l ) denote fitness values of the worst, the best, and the lth females, respectively. For minimum optimization problems, female spiders [like X l → ( t ) ] with low pheromone [ p h ( l ) ≤ 0.3 ] are difficult to mate with males, which will be replaced according to Eq. 5 to improve the quality of the population, as follows: X l → ( t + 1 ) = X b → ( t ) − 1 2 ( X r 1 → ( t ) − ( − 1 ) σ ⋅ X r 2 → ( t ) ) (5) where X r 1 → ( t ) and X r 2 → ( t ) are two different search agents randomly selected in the current population; σ stands for a random binary number, which is equal to either 0 or 1.

In order to acquire the best performance of frequency regulation, an optimal parameter extraction process of step start-up adaptive inertial droop controller based on BWOA is proposed in this section. On the one hand, this paper primarily concentrates on dynamic response characteristics of system frequency f AC after an increase in load illustrated in Figure 4, i.e., FFD, SFD, and TFD, that reflect local details, among which f FFD , f SFD , and f TFD represent the nadirs of FFD, SFD, and TFD, respectively. Specifically, cluster 1 is immediately put into operation after frequency event at t 0 (i.e., Δ t equal to 0 for WTs of cluster 1) while cluster 2 is launched when the rotor angular speeds of cluster 1 start to recover at t I , i.e., Δ t is 0 for WTs of cluster 2.

On the other hand, the integrals of frequency deviation and frequency variation rate ( S IoFD and S IoFVR , respectively) are also considered, which principally manifest global information of system frequency fluctuation, among which S IoFD and S IoFVR can be respectively expressed as S IoFD = ∫ | 50 − f A C | d t (6) S IoFVR = ∫ | d f AC d t | d t (7)

As indicated in Eq. 8, each indicator is assigned a coefficient according to the relative significance, and they are added up as the fitness function. F fit ( S IoFD ,   S IoFVR ,   f FFD ,   f SFD ,   f TFD ) = min { ( 10 A IoFD + 100 A IoFVR + 200 f FFD + 100 f SFD + 100 f TFD ) ⋅ 50 } (8)

Additionally, parameters of the controller for different WTs within the same cluster are identically designed for reducing the complexities of the problem and enhancing the search efficiency of BWOA. Therefore, optimized parameters can be further refined as X = [ k dr , i , 0 1 ,   k in , i , 0 1 ,   k dr 1 ,   k in 1 ,   k dr , i , 0 2 ,   k in , i , 0 2 ,   k dr 2 ,   k in 2 ] , in which the superscript 1 and 2 stand for cluster 1 and cluster 2, respectively. Besides, the lower and upper limits of optimized parameters are demonstrated in Table 1.

Finally, the parameter optimization pseudocode and flow chart of step start-up adaptive inertial droop controller with BWOA are explicitly exhibited in Table 2 and Figure 5, respectively, among which N and It _max denote the size of population and maximum iteration, respectively.

The dynamic frequency responses acquired by different methods under various load increase variations are elaborately depicted in Figure 6. One can easily observe that the system frequency based on BWOA optimization performs the least fluctuation in all load increase scenarios compared with those of other methods. Furthermore, BWOA-based parameter optimization method efficiently suppresses the FFD, SFD, and TFD of system frequency, which can dramatically ensure grid stability and reliability.

Furthermore, quantitative comparison results of frequency drop based on different strategies under three typical load increase situations are tabulated in Table 4, where the least frequency drop is highlighted in bold. In general, the trial-and-error method acquired the least FFD followed by BWOA and without FSC. However, the performances of the trial-and-error method on suppressing SFD and TFD were inferior to those of BWOA. In particular, compared with the trial-and-error method, the FFDs and TFDs obtained by BWOA are reduced by 10.9%, 9.8%, and 8.9%, as well as 12.8%, 12.2%, and 11% under 50 MW load increase, 150 MW load increase, and 250 MW load increase, respectively.

So as to confirm the effectiveness of control parameters acquired by BWOA in load decrease, Bus 1 of the test system experienced a 50-MW load decrease at 5 s, which was used as a simulation case. In this frequency event, cluster 1 (i.e., WT1, WT2, and WT2) is immediately started up at the fifth second while cluster 2 (i.e., WT4 and WT5) is launched at the 10th second. The simulation results are demonstrated in Figure 7 and Figure 8. It can be seen from Figure 7 that three control schemes caused 0.0641-, 0.0597-, and 0.0613-Hz frequency increases, respectively. Although frequency increase with the trial-and-error method is lower than that with BWOA at the sixth second, two obvious frequency fluctuations occur with the former one in the following procedures, which are caused by superabundant active power injection corresponding to Figure 8. Therefore, obtained optimal control parameters by BWOA are applicable to load decrease with low-frequency fluctuation.