An online identification method for establishing a microgrid equivalent model based on the hybrid particle swarm optimization butterfly algorithm

The frequency response model of a microgrid system is an indispensable tool for designing secondary frequency controllers and analyzing system frequency stability. Owing to the complexity of a microgrid system structure, time‐varying operation state, and difficulty in obtaining the inverter parameters, it is difficult to model the microgrid system. To address the abovementioned problems, this study constructs the equivalent model of a microgrid system based on mechanism analysis and then identifies the equivalent model parameters using an online identification method based on a hybrid particle swarm optimization butterfly algorithm according to the sampled operation data to obtain the specific equivalent model of a microgrid for the controller. This modeling method has strong adaptability and can yield accurate identification modeling when the operational structure of the microgrid changes.


| INTRODUCTION
Microgrid system modeling is an important method for studying the stability and optimal design of microgrid systems. 1 Modeling can be classified into mechanism and identification modeling. Mechanism modeling corresponds to building an accurate mathematical model of a microgrid system with known data regarding its topological structure and line parameters. [2][3][4] The main advantage of this method is its high accuracy, while its disadvantage is the complexity of the modeling process. Mechanism modeling is suitable for applications where the topological structure and line parameters are known; for example, the stability analysis of a single inverter. However, it is difficult to apply mechanism modeling when designing the parameters of the secondary frequency controller of a microgrid system. System identification modeling is an effective method for modeling external characteristics, and system modeling methods are classified into white-box modeling, 5 gray-box modeling, [6][7][8] and black-box modeling. 9 White-box modeling, also known as system parameter identification, corresponds to identifying system parameters when the operating mechanism of the system is completely known and is generally used in the design of controllers. Gray-box modeling is based on the premise that the structure and parameters of the system are partially known. With this framework, the model structure is generally obtained through theoretical analysis, which is then combined with experimental data to identify the system model. Black-box modeling is based on the premise of completely unknown system internal conditions. Accordingly, the model identification is based on the input and output data.
In the distributed generation field, the identification modeling approach was first applied to the modeling of photovoltaic systems and gradually extended to the identification modeling of distributed generation. [10][11][12][13] Moreover, some scholars have studied the identification modeling of distributed generation. For example, according to the characteristics of virtual synchronous generator machine control, Zaker et al. 14 and Yang et al. 15 proposed parameter identification methods based on the least-squares method and a genetic algorithm to identify the damping coefficient of a virtual synchronous generator machine. Luo et al. 16 identified the equivalent characteristics of grid-connected inverters; however, it focused on the parameter identification of a single inverter model and is not applied to the identification of the entire microgrid model. Yu et al. 17 proposed a black-box identification modeling method for a microgrid. By measuring the voltage frequency data at the PCC and total active and reactive power command values of the system, the voltage frequency response model of the microgrid can be identified using the least-squares method, which can effectively address the problem of difficulty in obtaining the parameters of microgrid equipment. However, identification in nonlinear black box systems contains a certain uncertainty. Therefore, it cannot be applied to the online optimization analysis of a microgrid system. Shi et al. 18 proposed an adaptive online identification method combining the Akaike information criterion using the ordinary least-squares method to update the system's optimal model order and parameters in real-time. However, this approach also suffered from problems of extensive computation, variable model order for each identification, and less apparent physical meaning.
In summary, for the entire microgrid, simpler modeling methods are required to build the microgrid model, and fast and accurate online identification methods are required to obtain the unknown parameters in the equivalent model in real-time. Therefore, this study proposes a method to construct the equivalent model of a microgrid system and identify the parameters directly online. First, the equivalent model of the grid is constructed through mechanism analysis. Then, the frequency and power data of the common connection points in the microgrid system are measured, and the hybrid particle swarm butterfly algorithm is applied to identify the equivalence parameters of the microgrid. Subsequently, the frequency response model of the microgrid is obtained.
2 | MICROGRID TOPOLOGY AND FREQUENCY CONTROL 2.1 | Introduction to typical microgrid topology A typical microgrid system topology is shown in Figure 1. The system includes multiple energy storage inverters, photovoltaic systems, and wind-power generation systems. In addition, it includes various loads and a microgrid central controller (MGCC). The intelligent circuit breaker in the microgrid and other devices with communication functions can obtain their operation statuses from the MGCC through the communication system and issue control commands to the inverters.
F I G U R E 1 Topological structure of a microgrid system. In a microgrid system, the inverters adopt the droop control framework to achieve frequency stability. In the droop control framework, the relationship between the frequency and active power is stated as follows: Where P Oi is the set power base point value, f Oi is the set power base point frequency, P is the active power of the inverter, K m is the active power frequency droop coefficient, and f droop is the reference frequency.

| Microgrid secondary frequency control
Because droop control is a frequency difference control approach, the frequency in the microgrid varies with respect to the load. PQ control is adopted for photovoltaic power generation and wind power generation in microgrid systems, which do not participate in frequency control. The current response speed of photovoltaic and wind power generation is also very fast, which has no impact on the secondary frequency control of the system. The impact of power fluctuation on the system frequency after connecting to the microgrid is similar to the impact of load switching on the system frequency. The secondary frequency control strategy of the power system and the frequency control block diagram of the microgrid system are shown in Figure 2.
First, the MGCC reads the frequency signal from the point of common coupling (PCC), calculates the difference between the frequency feedback and given reference, and calculates the total power adjustment value through the controller; then distributes it according to the distribution coefficient K ai of each inverter. The inverter can adjust the droop curve according to the adjusted power based on the power and output power difference generated by the MGCC. In Figure 2, K mi denotes the droop coefficient of the ith inverter, K zi denotes the droop output power characteristic parameter of the ith inverter, being K πU E Z = 2 / zi 1 1 , P Oi denotes the reference power of the ith inverter, f Oi denotes the reference frequency of the ith inverter, and τ is the lowpass filter time constant of the inverter output power. To ensure the response speed, the sampling period of the MGCC is usually set to tens of milliseconds. The delay link can be approximated as an inertial link through Taylor expansion. Therefore, τ corresponds to the communication delay of the MGCC sampling frequency. To satisfy the requirements of power-sharing, the power distribution coefficient is usually selected according to the following equation: The microgrid system control block diagram shown in Figure 2 can be transformed into the simplified structure shown in Figure 3. As shown in Figure 3, the transfer function between the total output power P O of the inverter and the total regulated power ΔP O of the system is very complex. It is necessary to further create a simplified equivalent model of the microgrid secondary frequency control system. K mx indicates the equivalent droop coefficient of multiple inverters participating in the frequency control, and K zx indicates the output power characteristic parameters. During the actual modeling process, it is not necessary to accurately calculate the equivalent parameters K mx and K zx but only to approximate them through circuit analysis. According to the droop characteristic equation, the equivalent droop coefficient K m of the inverter can be expressed as follows: Through circuit analysis, the slight difference in the output voltage caused by the line impedance is ignored. Subsequently, the output power characteristic parameters K z are stated as follows: Therefore, the droop parameters and output power characteristic parameters of the equivalent single inverter are equal to their approximate values K m and K z , respectively, The microgrid secondary frequency equivalent model shown in Figure 4 was obtained after simplification and substitution.
To identify the parameters of the equivalent model, gray-box identification was conducted directly using the collected active power input of the microgrid and frequency data. We can take P O of Figure 4 as input and f PCC as output, thus obtaining Figure 5, so that we can obtain Equation (7).

| Equivalent parameter identification of the microgrid
According to the analysis presented above, it is necessary to identify the equivalent droop control parameter K m of the microgrid and select the appropriate sampling step to collect the frequency, inverter output voltage, and current. First, the equivalent transfer function (7) is expressed as the following differential equation: The recursive equation in the discrete domain can be obtained using the three-point interpolation derivation equation as follows: Where i indicates the ith data of the variable in the discrete domain, h is the sampling step, K m is the equivalent coefficient, and τ is the sampling period that corresponds to 20 ms. K p and K i are the proportional and integral controller gains of the system without adaptive optimization during sampling. The mean-square error method was used to measure the deviation between the simulated and true values; the results were used as the assessment standard to evaluate the parameter identification results. When the number of sampling points is n, the mean square error of the output frequency of the microgrid equivalent model with multiple inverters in parallel is stated as follows: Where f ′(i) represents the frequency value calculated using the identified equivalent coefficient. f (i) represents the frequency value sampled at the PCC. Fitness takes the sum of the difference between f (i) and f ′(i). When the sum is smaller, it indicates that the calculated frequency value is closer to the actual output value, which means that the identification result becomes more accurate.

| PSO algorithm
PSO is an evolutionary computation technique. 19 It originated from studying the behavior of flocks of birds foraging for food. The basic idea of the PSO algorithm is to determine the optimal solution through collaboration and information sharing between individuals in a population. Particles have two properties: velocity, which represents the speed of movement, and position, which indicates the direction of motion. The initial position and velocity of each particle in the search space are randomly initialized. The position and velocity of a particle are updated as follows: The v i t and v i t+1 indicate the velocity of the ith particle at t and t + 1 iterations, respectively. p best denotes the optimal position experienced by the particle itself, while g best denotes the optimal position experienced by the entire particle population. At each iteration, the adaptation value of the new position assumed by each particle is compared with p best and g best . In general, c 1 = c 2 = 2 and rand 1 and rand 2 are random numbers in the range (0,1). The two formulas presented above are used as the basis for the standard form of the PSO.

| Butterfly optimization algorithm (BOA)
Arora proposed a BOA, which is a novel metaheuristic intelligent optimization algorithm. 20 The algorithm was inspired by the foraging and courtship behaviors of butterflies. Butterflies acquire and analyze the scents in the air to determine the possible direction of their target. In this framework, a butterfly perceives odors with the following intensity: In the equation, f ra denotes a butterfly's perception of fragrance, c indicates the sensory morphology coefficient, I indicates the stimulus intensity, which corresponds to the function fitness value F(i) in Equation (10), α indicates the power index based on the degree of aroma absorption, within a range of [0,1].
BOA determines the global and local search results of the algorithm according to the switching probability p. The location update equation is stated as follows: YONG ET AL.

| 1623
Where x i t+1 indicates the spatial position of the ith butterfly at moment t, g* is the best position of all butterflies in the current iteration, and x j t and x k t indicate the spatial position of the jth and kth butterflies at the moment t, respectively. The value of r is a random number between [0,1] and f i is the fitness value of the ith butterfly.

| Butterfly colony initialization
When identifying the parameters of the microgrid equivalent model, aiming to address the problems of low convergence accuracy and ease of falling into local optimization experienced by the basic butterfly algorithm, a hybrid particle swarm optimization butterfly algorithm (HPSOBOA) is proposed to identify the microgrid equivalent model parameters.
Let the equation of the randomly generated initial solution in the d-dimensional search space be Here, X i is the spatial position of the ith butterfly (i = 1, 2, 3, …, N) in the butterfly population, N indicates the number of initial solutions, L b and U b indicate the upper and lower bounds of the search space, respectively, and o indicates the random number matrix between (0,1).

| Global search for algorithms
The activities of butterflies looking for food and mating partners in nature can be performed both globally and locally. Therefore, a switching probability p is set to convert the ordinary global search to a dense local search. In each iteration, HPSOBOA randomly generates a number r within [0,1] and compares it with the switching probability p to decide whether to conduct a global or local search. When r < p, a global search is initiated, which iterates according to the first line in Equation (16) as Where ω denotes the inertia weighting factor and X i t indicates the spatial position of the ith butterfly at moment t. Let V i t+1 and V i t indicate the velocity of the ith particle at moment t + 1 and t, respectively. Where c 1 = c 2 = 0.5, ω denotes the inertia weighting factor, and r 1 and r 2 are random numbers in (0,1) Update the spatial position of the particle at the moment t + 1 based on X i t and V i t+1 , where V i t+1 denotes the velocity of the ith particle t + 1, which can be calculated using Equation (17).

| Local search of the algorithm
When r > p, a local search is initiated, which iterates according to the second line in Equation (19) as Where X i k and X j t−1 denote the jth and kth butterflies randomly selected in the solution space and ω denotes the inertia weighting factor. Let Update the spatial position of the particle at the moment t + 1 based on X i t and V i t+1 .

| Control parameters ω and c
The meanings of ω and c mentioned above are as follows.
The inertia weight factor ω directly impacts the flight speed of the HPSOBOA algorithm, which can be used to adjust the global and local search capabilities of the algorithm. 21 The adaptive adjustment strategy used in this study is stated as follows: , and T max is the maximum number of iterations. 20 Theoretically, the sensory morphology coefficient c can take any value in the range [0,∞]. However, its value is determined by the optimization problem specificity during the butterfly algorithm iteration. 20 The sensory modalities occurring in the optimal search phase of the algorithm are expressed as follows: Here, T max is the maximum number of iterations in the algorithm, and the initial value of parameter c t is set to 0.01. 19

| Steps of the parameter identification algorithm
The steps of the parameter identification algorithm used in the microgrid equivalent model based on the HPSOBOA are as Figure 6. a. Initialization: Set the parameter values required by the algorithm, initialize the parameter K mi to be identified, and randomly generate a matrix of N × b according to the upper and lower limits of the search space, where N is the number of individuals and b is the individual dimension. b. Calculate the fitness value and determine the optimal flower nectar source position of the identification parameter butterfly K mi . Calculate the fitness value of the parameter butterfly K mi to be identified according to the objective function based on the error sum. The position of the minimum parameter butterfly K mi to be identified corresponds to the flower nectar source position. Then, calculate the fragrance size of the parameter butterfly K mi to be identified according to Equation (13). Because the fitness function is used to find the minimum value, the stimulation intensity I takes the reciprocal 1/F of the objective function. Accordingly, the more accurately the parameter is identified, the greater the fragrance of butterfly K mi . c. Update the butterfly K mi position of the parameter to be identified: Judge whether the current iteration performs a global or local search according to the switching probability p, and update and determine the position of each butterfly. d. Repeat the steps to calculate the fitness value and update the position of the parameter butterfly to be identified until the predetermined iteration end condition is reached. Subsequently, terminate the algorithm and output the optimal solution.

| SIMULATION EXPERIMENTS
In this study, a series of parallel-output microgrid inverters with a rated power of 100 kW is used to verify the accuracy of the proposed parameter identification method. The three-phase load parameter is set to 30 kW. All three inverters operate in the P-f droop mode with droop coefficients of K m1 = 8 × 10 −6 , K m2 = 5 × 10 −6 , and K m3 = 3 × 10 −6 . The model is provided in Figure 7.

| Comparison experiments conducted on the identification results of multiple algorithms
To verify the robustness and effectiveness of the proposed algorithm, 20 independent comparative experiments of model parameter identification were carried out using the abovementioned microgrid equivalent model and HPSOBOA, PSO, and BOA algorithms. The maximum number of iterations in the parameter identification experiment was set to 500; the number of individuals in the population was set to 30. The parameter settings for each algorithm are presented in Table 1.
Calculate the average fitness convergence value of 20 independent experiments: Where F t best denotes the average fitness convergence value, reflecting the average change law of the global optimal fitness value in the iterative process, M denotes the number of experiments, and f t best denotes the global optimal fitness value at moment t for each independent experiment. The variation curve of the root-mean-square error objective function during the iterative process of the three algorithms is shown in the following Figure 8.
Observe from the figure that HPSOBOA has a smaller fitness value with the same iterations and is more likely to achieve the accuracy required for fitness. Similarly, with the same fitness value, HPSOBOA requires fewer iterations. The HPSOBOA algorithm used to identify the parameters of the microgrid equivalent model can meet the accuracy requirements faster compared with the other algorithms, while achieving a higher convergence accuracy. This indicates that the HPSOBOA algorithm attains a fast convergence speed, can reduce the number of iterations, and has better stability in high-numericaloptimization problems. Because the time of the microgrid system simulation is much longer than the algorithm runtime, reducing the number of iterations is crucial in reducing the total time required for parameter identification.

| Equivalent parameter identification experiment
To verify the ability to identify the equivalent coefficients of the HPSOBOA algorithm proposed in this study, the following simulation scheme was designed and compared with the ordinary least squares used to identify the parameters of microgrids. The simulation droop parameters of the three inverters were set to K m1 = 8 × 10 −6 , K m2 = 5 × 10 −6 , and K m3 = 3 × 10 −6 . Substituting these values into Equation (4), the theoretical value is obtained as K m = 1.52 × 10 −6 . A power disturbance is applied, with sampling step h of 10 −3 s. The simulated sampling data are provided to the least-squares and HPSOBOA programs. The parameter settings of HPSOBOA are the same as those in described Section 4.1. To verify the accuracy of the equivalent model identified based on the PSO hybrid butterfly algorithm, a simulation scheme of the microgrid under the frequent equipment switching condition is designed, where "1" indicates that the inverter is connected to the microgrid and "0" indicates that the inverter is disconnected from the microgrid. The settings are listed in Table 2.
For example, when the unit is switched on and off, the MGCC starts the HPSOBOA identification program to obtain the real-time equivalent droop coefficient of the microgrid. During the identification process, optimization is conducted based on the fitness equation. The specific identification results are as follows.
The equivalent coefficient is stated on these three time stamps. The unit is switched on and off at 1 s and 3 s. The results are shown in Figure 9. Abbreviations: BOA, butterfly optimization algorithm; HPSOBOA, hybrid particle swarm optimization butterfly algorithm; PSO, particle swarm optimization.
The system structure is the same between 2 and 4 s. The first converter operates alone for 4 and 6 s. The remaining two converters are connected in parallel. The corresponding equivalent parameter identification results are the same as well. The identification results obtained using the traditional least-squares method demonstrate that the difference between the identification results and HPSOBOA is small.
The identified equivalent parameters were substituted into the identification model. The active power input P O was applied to the model; accordingly, the output f PCC of the identification model and measured production were compared with both algorithms. The results are shown in Figure 10.
Compared with the traditional least-squares method, the model frequency output identified by the HPSOBOA algorithm agrees well with the actual output frequency waveform of the grid. This reflects the actual frequency output of the microgrid.

| Advantages of the dynamic experiment
To further verify the adaptability of the improved butterfly algorithm proposed in this study, namely, the online identification effect when the microgrid structure changes. When the microgrid structure varies, the droop coefficients of the three inverters are K m1 = 8 × 10 −6 , K m2 = 5 × 10 −6 , and K m3 = 3 × 10 −6 . The load of the microgrid and the number of inverters connected are changed, and the identified equivalent parameter values F I G U R E 8 Comparison between the optimization abilities of the three algorithms. F I G U R E 9 Hybrid particle swarm optimization butterfly algorithm (HPSOBOA) and ordinary least squares (OLS) identification results.
F I G U R E 10 Transfer function model fitting curves obtained from the active power step measured frequency simulation. HPSOBOA, hybrid particle swarm optimization butterfly algorithm; OLS, ordinary least squares.  Table 3.
The waveforms of the entire experimental process are shown in Figure 11. The butterfly algorithm of hybrid PSO is used for online identification; subsequently, the output frequency of the identification model was obtained. To test the adaptability and accuracy of the identification model, its output frequency waveform was compared with the measured output frequency waveform of the microgrid. A comparison between the measured output and model output is shown in Figure 11.
Observe that the output frequency amplitude of the simulation model can accurately reflect the real output of the microgrid. Especially when the frequency changes, it can be accurately tracked, which underlines that the identification algorithm can adapt well to changes in the microgrid structure.

| CONCLUSION
In this study, the butterfly algorithm of hybrid PSO was used to establish the frequency response gray-box model of the microgrid. This study focuses on the parameter identification algorithm and method of establishing a microgrid gray-box model using the corresponding transfer function. Compared with traditional microgrid parameter identification, the proposed method improves the accuracy and identification speed. Finally, by comparing the model output with the measured output, the identification model can quickly and accurately reflect the frequency characteristics of the microgrid by relying on the HPSOBOA algorithm when the microgrid structure changes. This indicates that the gray-box identification modeling method proposed in this study has wide adaptability to realize an online system identification.