Power Grid Frequency Regulation Strategy for Photovoltaic Plant Based on Multi-Objective Harris Hawks Optimization

To settle the issue of balance between two objectives, i.e., photovoltaic (PV) power station output power maximization and frequency regulation (FR) signals response, a novel PV reconfiguration strategy is proposed in this work, which maximizes the output power through PV reconfiguration, and meanwhile utilizes the energy storage system (ESS) to decrease the PV plant generated power’ deviation from FR signals. Above all, a model of PV-storage power station reconfiguration is designed to minimize the power bias of both rated power and FR signals. Then, the multi-objective Harris hawks optimization (MHHO) is used to obtain the Pareto front which can optimize the above two objectives due to its high optimization efficiency and speed. Subsequently, the optimal compromise solution is selected by the decision-making method of VIseKriterijumska Optimizacija I Kompromisno Resenje (VIKOR). Aiming to substantiate the efficacy of the proposed technique, the case studies are carried out under partial shading condition (PSC) with constant and time-varying FR signals. The simulation results show that, compared with the situation without optimization, the power deviations of the two objectives are reduced by 25.11 and 75.76% under constant FR signals and 23.27 and 55.81% under time-varying FR signals by proposed method, respectively.


INTRODUCTION
Nowadays, major countries in the world are committed to developing the application technology of renewable energy resources to tackle the incoming traditional energy crisis. With respect to its abundant reserves and environmentally-friendly quality, solar energy has become the most widespread energy resource (Gan et al., 2019;Mahidin et al., 2021). In general, photovoltaic (PV) arrays operate at the maximum power point by maximum power point tracking (MPPT) (Salam et al., 2013;Yang et al., 2016) technique to maximize the outputs of PV station. However, in the field application of large-scale PV power generation system, the bypass diodes of PV panel will be activated under partial shading condition (PSC), resulting in inconformity of output characteristic between various PV array rows (Sai Krishna and Moger, 2019;Zhang et al., 2021a). And further consequence of multi-peak with power-voltage (P-V) curve is caused which can prominently decrease PV plant power generation and make it more difficult for MPPT. Besides, PSC could lead to hot spot phenomenon which may break PV panels and even cause a fire. In addition, PSC is mostly attributed to clouds, building shadows, dust, birds' feces, manufacturing defect and non-uniform aging of PV components (Yang et al., 2018a;Laudani et al., 2018).
For the purpose of mitigating the effects of PSC as well as enhancing the power generation of PV station, altering the topological structure of PV components is employed, which could equivalently disperse the originally-centralized shadow on the array surface to the whole array. In recent years, researchers have developed multiple methods of PV arrays reconfiguration, which are classified mainly into static reconfiguration techniques and dynamic techniques (Kaushika and Gautam, 2003;Yang et al., 2021).
Three different topological structures of PV arrays, i.e., series-parallel (SP), bridge linked (BL), and total-cross-tied (TCT), are studied and compared under PSC in literature (Kaushika and Gautam, 2003;Ajmal et al., 2020), of which the performance of TCT is ascertained the best topological configuration theoretically for its robustness. But TCT has higher wiring losses for its longer cable. Literature (Rani et al., 2013) reconfigures the PV components based on Sudoku patterns, and a variant of sudoku-based arrangement (Horoufiany and Ghandehari, 2018) is proposed to tackle with mutual shadows (Horoufiany and Ghandehari, 2017;Yang et al., 2020a). By introducing mixed integer quadratic programming, literature (Shams el-Dein et al., 2013) utilizes the mathematical formulation to dissipate the shadow of PV array. TomTom puzzle divides the PV array into irregular pieces with particular calculation rule assigned for each piece, and then accomplishes the reconfiguration of PV array (Tatabhatla et al., 2019;Yang et al., 2019). The essential characteristic of static reconfiguration is that only the physical location of the PV module but not the electrical connection will change. Moreover, as a one-time reconstruction technique, static method is not able to automatically alter reconfiguration schema to meet various shadow pattern. Therefore, the PV system with static method must artificially constructs a new configuration when a new shadow appears. In practice, shadow patterns frequently change, hence static method is not universally applicable, nor is it suitable for real-time scenarios.
Dynamic PV array reconfiguration techniques come into being in response to the aforementioned drawback of static techniques. Different from static methods, dynamic methods can dynamically change the electrical connection of PV modules in real time by switching matrix according to the shadow change (La Manna et al., 2014;Yang et al., 2018b). Literature (Vicente et al., 2015) applies rough set theory to PV array reconfiguration. Literature (Sanseverino et al., 2015) builds a control system embedded with Munkres algorithm to alter the switches layout. Feed forward artificial neural network is utilized to alter the switch matrix by training irradiance data and shortcircuit currents (Karakose et al., 2014). Nowadays, meta-heuristic algorithms has become a prevailing tendency in applicating to PV array reconfiguration, such as genetic algorithm (GA) (Deshkar et al., 2015), flow regime algorithm (FRA) , particle swarm optimization (PSO) (Babu et al., 2018), grasshopper optimization algorithm (GOA) (Fathy, 2018), butterfly optimization algorithm (BOA) (Fathy, 2020), artificial ecosystem-based optimization (AEO)  and so forth.
However, existing researches take little account of the combination of PV array reconfiguration and power grid frequency regulation (FR) (Xi et al., 2018), nor do they take the power fluctuation balance of energy storage system (ESS) (Jin et al., 2017;Zhou et al., 2021) into consideration. Hence, this work proposes a novel technique of PV array associating power grid FR, to promote the dynamic-response performance of the entire regional power grid.
Given the confliction between power outputs maximization and minimizing FR signals' deviation from the output power, this work employs the efficient multi-objective Harris hawks optimization (MHHO) to obtain the Pareto front promptly, then Visekriterijumska Optimizacija I Kompromisno Resenje (VIKOR) strategy is adopted to determine the optimum compromising solution based on aforementioned Pareto front.

Modeling of PV Array
Each PV module is formed by multiple of PV cell strings connected in parallel, while PV array is comprised of several PV modules connected in series and parallel. Accordingly, multiple PV arrays make up a PV system. The most frequently employed model is the single diode PV cell model (Jacob et al., 2015;Humada et al., 2016;Yang et al., 2020b), illustrated as Figure 1A. Calculated by Kirchhoff's Current Law (KCL), the operation current of PV cell I cell is given below: where I Light denotes the photo-generated current; I d is the current of diode; I sh represents the current of shunt resistor. By replacing the I d and I sh , the Eq. 1 can be rewritten as Eq. (2).
where I o and q denote the leakage current and electronic charge of the diode, respectively; T c , V cell , R sh , R s , σ and B are the operation temperature and voltage, shunt and series resistance of PV cells, ideality factor and Boltzmann's constant.
The output current I m of one PV module comprised of n s PV cells in series can be calculated as follows: where I m denotes the PV module current output; V m represents the operation voltage of PV module; R SH and R S separately refer to the parallel and series resistor of PV module. Moreover, I Light represents the photo-generated current of PV module which is caused by photoelectric effect, and its formula is given below: where G and G 0 , I stc , T c and T 0 , K sc represent the actual and reference irradiance, the photo-generated current under standard environment condition, the operating and standard temperature, and the short-circuit current factor, respectively. Thus, as demonstrated in Figure 1B, it can be further deduced the output current I a of the PV component formed by N s ×N p PV modules can be described as: where V a and I a respectively refer to the operation current and operation voltage of PV components.

Modeling of PV Array in TCT Structure
Because of the better stability of TCT topology under PSC, it is widely applied to the theoretical research of PV reconfiguration (Bingöl and Özkaya, 2018). As depicted in Figure 2, the 9 × 9 PV array is connected in TCT configuration, the output voltage and current of which can be calculated by: where V S represents the entire voltage of PV system; V ai represents the row voltage at the i th row; I S represents the entire current; and I ij is the flowing current at the i th row and the j th column.

Optimization Objective
This work features two sets of assignments: one is to reduce the power deviation y 1 of real PV power generation average value from the rated PV power; while the other aims to minimize average value y 2 of the power difference between the outputs of PV station and FR signals. Meanwhile, the peak-to-valley deviation of the power curve is taken into consideration to control its fluctuation. In consequence the proposed objective function can be computed as follows: where P rate pv represents the rated output power of PV plant; T denotes the duration of FR; P pv (t), P ESS (t), P FR (t) respectively denote the generated power of the PV station with ESS, the charging power of ESS, and the FR signals in the t th control time period.

Constraint Condition
Each PV module can only switch to other row in the same column to reconfigure PV system via switching matrix, which is called constraint condition of the electrical switch states. The power and capability constraints should be satisfied in the ESS, shown as follows: where P max ESS , P min ESS , SOC max ESS , SOC min ESS severally represent the maximum and minimum charging power, and the maximum, minimum state of charge (SOC) of ESS; η ch and η dis are the charging and discharging efficiency; Δt denotes the control period; and E Ess denotes the rated capacity of ESS.

OPTIMIZATION FOR PV-STORAGE STATION RECONFIGURATION BASED ON MHHO AND VIKOR Design of MHHO
MHHO simulates the cooperation mechanism and pounce style in the hunting of Harris hawk's predation. The purpose of this research is to integrate the non-dominated rank sorting with MHHO, evolving into a multi-objective optimization algorithm.

Exploration Stage of MHHO
In MHHO, the Harris hawks can be seen as the candidate solution, and in each iteration process the optimal solution is considered as the prey. Harris hawks apply two types of inhabitation tactics with equally mathematical probability, the specific model of which is demonstrated as follows: where X(t + 1) and X(t) respectively denote the hawk's position vector in the next iteration and in the present; X prey (t) is the position vector of the prey; r 1 , r 2 , r 3 , q are random numbers between 0 and 1; U b and L b represent search bounds; X rand (t) denotes the stochastic individual in the present population; and X m (t) denotes the population's average position vector.

Exploitation Stage of MHHO
The flock of hawks round up the prey softly or strongly from all directions by tracing the left escape energy of prey. MHHO algorithm introduces the parameter of E: the flock performs soft besiegement based on Eq. 13, 14 if |E| < 0.5, otherwise the hard encirclement would be executed based on Eq. 15.
where J 2(1 − r 5 ), and J represents the prey's leaping distance when they act on escaping; and r 5 is a digit randomly generated by the computer between 0 and 1.

Multi-Objective Non-dominated Sorting and Archive Management
MHHO selects an optimal Pareto front via non-dominated rank sorting, then introduces archive management to control Pareto solutions. When a new solution is created, it is to be compared with the ones in the archive. Were it dominated by those, the new solution would not be collected into archive; Were a few individuals in the archive dominated by new solution, those would be replaced by the new one; if it is none of the above, this new solution would also be collected. Additionally, when the quantity of solutions exceeds the archive' capacity, the roulette shall roll out the more crowded solutions.

Fitness Function
On the basis of the penalty function technique, the fitness function of MHHO is calculated by Eqs 16, 17.  where Pe(t) denotes the penalty part of the restraint conditions in Eq. 10 in the t th dispatching period of FR signals; and ϕ denotes the penalty factor that usually employs a large positive number.

Design of VIKOR
The optimal compromise solution can be impersonally chosen by VIKOR (Lin et al., 2021) decision-making approach from the Pareto front solutions obtained by MHHO, the design process of which is elaborated below: Step 1: Define the plus-minus ideal solutions on the basis of the obtained Pareto front, as given: where v i j is the corresponding i th Pareto numerical solution of the j th objective function; and v + and v − respectively represent the plus and minus perfect solutions.
Step 2: Calculate the values of the cluster effect and single regret: where φ j is the weighting factor of the j th objective function; and U i , R i denote the cluster effect value and single regret value corresponding to i th Pareto front solution, separately.
Step 3: Determine the integrated valuation of all the Pareto front solutions: where c means the weighting factor of the cluster effect value; Q i denotes the integrated valuation value of the i th Pareto front solution; U + , U − denote the minimum and maximum values of cluster effect; and R + , R − represent the minimum and maximum values of single regret.
Step 4: Find the optimal eclectic solution on the basis of integrated valuation value, as follows: x opti arg min i 1,2,...,n where x p i represents the i th Pareto solution; x opti is the optimal eclectic solution.

Procedure Steps
With regard to the PV station array reconfiguration, the procedure that MHHO executes is demonstrated in the following Table 1.

CASE STUDIES
The simulation experiment conducted herein designs a 40 MW PV power plant comprised of 50 identical 9 × 9 subsystems, to test the MHHO performance in the experiment of PV power plant' response to FR signals. The parameters of the subsystems are derived from literature (Zhang et al., 2021b), and the specific parameter settings of PV module are tabulated in Table 2. This work supposes all PV arrays' operation temperature to be at 25°C, and the shadow distribution of each sub-system to be consistent. Figure 3A illustrates sub-system's shadow distribution in 10 min, while Table 3 gives the main parameters of ESS.
The number of agents (N pop ) is 200. Meanwhile the termination number of iterations (Max itera ) is set as 100. Moreover, all of the weighting factor in VIKOR are set as: φ j 0.5. All simulation tests are conducted on Matlab R2019b.

Constant FR Signals
To testify the effectiveness of the proposed algorithm under the constant FR signals, the constant FR signals are designed as follows: The shadow distribution of PV array after reconfiguration under constant FR signals is depicted in Figure 3B. Obviously, MHHO commendably disperse the originally-centralized shadow to the whole array. Meanwhile, the optimal compromise and eclectic solution obtained by VIKOR exhibits more balance among two objectives than other Pareto front solutions, as it is shown in Figure 4A. In particular, the results of optimal eclectic solution are: the power deviation between the average generated power of PV array and the rated value is 9.128MW; the average power deviation between the output power of PV-storage station and constant FR signals is 2.903MW. Figure 4B illustrates the output power with optimization under constant FR signals of PV array, PV power plant, and ESS. With the help of ESS, the output power of PV plant is higher than that without optimization and gets closer to FR signal. The power outputs of PV power plant before and after optimization are given in Figure 4C and Figure 4D. After optimization by the proposed method, the objective f 1 drops sharply by 25.11%, while the objective f 2 drops by 75.76%. In addition, the output power with optimization experiences less power fluctuation.

Time-Varying FR Signals
The time-varying FR signals are designed to testify the effect of the proposed approach under the condition of it, as follows: (27) Figure 3C shows the irradiation distribution of PV array after reconfiguration under time-varying FR signals, of which the shadow distribution is more uniform compared to that without optimization. In the meantime, the compromise solution obtained by VIKOR shows better equilibrium among multiple targets than other Pareto front solutions, as demonstrated in Figure 5A. Particularly, the results of optimal compromise solution are: the power deviation between the average generated power of PV array and the rated value is December 2021 | Volume 9 | Article 806263 9.353MW; the average power deviation between the output power of PV-storage station and constant FR signals is 4.708MW. Figure 5B depicts the output power with optimization of PV array, PV power plant, and ESS. It can be deduced that due to the coordination of ESS, the output power of PV plant gets closer to FR signal. Figures 5C, D gives the results of PV power plant's power output before and after optimization. After optimization by the proposed method, the general deviation between the output power of PV plant and FR signals drops sharply by 55.81%, while the power deviation between the PV plant output power and rated power drops by 23.27%. Compared with constant FR signals, the optimization effect of objective f 2 slightly drops while that of objective f 1 is close. And the possible reason is that time-varying FR signals are more complicated. Also, the output power with optimization experiences less power fluctuation.

CONCLUSIONS AND PERSPECTIVES
Founded on PV array reconfiguration, a new technique of PV station participating grid FR is presented in this work. The proposed technique contributes to effective power response of PV plant to FR signals, and is consequently conducive to improving the frequency stability of the regional grid. The main contributions of this work are summed up as follows: 1) Not only can the proposed MHHO algorithm optimize two conflicted objectives simultaneously, but reconfigure PV array for FR expeditiously. Moreover, MHHO can find a series of solutions that fit simultaneously two objectives to meet the different dispatch needs; 2) To entail the balance of each objective function, VIKOR decision-making technique is adopted to select the optimal eclectic solution, which gives a guarantee of the balance between PV power station's economical operation and its security; 3) By virtue of the proposed method of PV array reconfiguration coordinated with ESS, the PV station is able to track the maximum power point as possible under any irradiance distribution, reducing the abandoned solar energy. In the meanwhile, the proposed method increases both economic benefit and security of the PV station.
However, it needs many independent runs to find the optimal solution due to the high randomness and local convergence of MHHO. Moreover, the implementation of ESS increase the cost of FR response. In addition, the proposed strategy cannot work in the condition of no shadow. Besides, The future work will focus on the following perspectives: • Except for simulation test, hardware-in-the-loop (HIL) experiments should be performed to examine the practical effect of proposed method; • More studies on larger scales of PV array are urgent to testify the universal applicability of proposed approach. • Weighting parameters are difficult to adjust which significantly affect the convergence speed of metaheuristic algorithms and the spatial distribution of the optimal solution and hence deserve deeper research; • In addition to meta-heuristic algorithms, hybrid algorithm and AI algorithm should be applied in this field of PV array reconfiguration participating grid FR; • In addition to VIKOR, more attention should be paid to other decision-making technique to find more reasonable compromise solution.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.