Harris hawks optimization algorithm for load frequency control of isolated multi-source power generating systems

2.1 System overview

In this work, LFC is examined in an isolated multi-source power system. The system under consideration consists of a reheated thermal-turbine power source, a hydro-turbine power source and a gas-turbine power source [5]. Figure 1 depicts the equivalent small signal perturbation block diagram. Under normal conditions, there is no inequality between production and demand and the following relationship holds [23]:

where the power of the steam generation source P _Gs = K _s P _G , power of the hydro generation source P _Gh = K _h P _G and power of the gas generation source P _Gg = K _g P _G . The K _s , K _h and K _g are the sources’ participation factors (spf), which may be used for systems with more than one type of source per control area, like the one examined here. In some studies, source participation factors divide the power of each area’s loads equally, so that all units serve the corresponding portion of the demand, but this is not true in real systems whereas the demand is divided according to the availability of the generators comprising the source [5]. However, in any case, the following constraint applies [23]:

Figure 1:

Isolated multi-source power system under study (uncontrolled).

In this particular system, the total power is 1750 MW, of which the steam, hydro and gas sources contribute 1000 MW, 500 MW and 250 MW respectively (ratio of 4:2:1). In this loading condition, the participation coefficients are derived and given in the Appendix along with all the other system parameters. The power system subsystem in Figure 1 can be modeled under the assumption that the power of the turbines is balanced by the power of the synchronous generators. The demanded load power is subject to continuous variations ΔP _L (s), which must be covered by the generator’s power. The difference between generated power and demand creates a variation in the frequency Δf(s) of the system. The equilibrium equation is then given by the following relationship:

This equation can be written as [23]:

where K _ps is the power system gain constant (in Hz/pu MW) and T _ps is the power system time constant which can be obtained by [5]:

where H is the system inertia constant, D is the damping factor and f ₀ is the system rated frequency.

2.2 System inherent non-linearities

In real power plants, restrictions apply (in the form of non-linearities) known as generation rate limits which exist on the rate of change in power generation and also in the governor’s action. The first is called generation rate constraint (GRC), while the latter is called governor’s dead band (GDB). In literature, these two restrictions are usually met in steam and hydro power sources. Also, in a steam source, boiler dynamics (BD) can be considered in LFC studies, which also impose increased non-linearity on the system. The block diagrams of these constraints are shown in Figure 2 and are explained next.

Figure 2:

With respect to Figure 1: (a) boiler dynamics block diagram (BD), (b) governor dead-band dynamics (GDB), and (c) generation rate constraint block diagram (GRC).

2.3 Control strategy and optimization function

The main step in testing the controller design with an optimization algorithm is to choose appropriate objective functions. The integral absolute error (IAE), integral-squared error (ISE), integral-time-squared error (ITSE), and integral time absolute error (ITAE) are examples of integral-based objective functions. In contrast to the ITAE, the ISE, IAE and ITSE criteria require longer settling times and cannot always perform within a desirable stability margin [25, 26]. The literature has shown that the ITAE is more effective than other performance criteria. The ITAE is used as the optimization objective function for the controller gains in this study, and it is given by:

where T _sim is the time range of the simulation and Δf is the frequency deviation of the system.

As mentioned in the introduction Section, there are many configurations of controller types that may be used in LFC studies. In the current study, the system has been modeled using a PI, PID and FOPID controller for each plant in MATLAB/Simulink environment.

3.1 Exploration phase

The candidate solutions in this algorithm are hawks, and at each phase the ideal candidate solution is regarded as the desired prey. The algorithm predicts that hawks will randomly roost in some places and wait until they spot a prey using two different tactics. Assuming equal probability p between (0, 1) for each tactic, the hawks either roost according to where the other group members are located (so that they are close to them when they attack) and the position of the hare, a first strategy modeled for p < 0.5, or they roost randomly in tall trees (random locations within the group area), a second strategy modeled for p ≥ 0.5.

where v ₁, v ₂, v ₃, and v ₄ are random integers between (0, 1) that are used in every iteration, the hawks’ position vector at the current and the next iteration is denoted by Y(k) and Y(k + 1), Y _prey(k) is the position of the hare, Y _rnd(k) is a hawk (randomly chosen) out of the current population and Y _a represents the mean location of the current hawk population. The lower and upper bounds are denoted by L _b and U _b , respectively [28]. The hawks’ average location is determined by the following equation:

where M represents the total population size and Y _i (k) indicates each hawk’s location at iteration k [29].

3.2 Exploration to exploitation transition

Depending on how the prey flees, the HHO algorithm may switch between various exploitation behaviors and exploration ones. The following equation models how prey energy is drastically reduced during escape behavior:

where P denotes the escape prey’s energy, K refers to the maximum iteration number and P ₀ denotes the initial state of its energy inside the interval [−1, 1]. When the value of P ₀ decreases from 0 to −1, the hare loses energy. When the P ₀ value increases from 0 up to 1, the hare is getting stronger. The escape potential energy P during iterations has a declining trend. When |P| > 1, the hawks look for other areas to explore locations with hares, and the algorithm executes the phase of exploration. When |P| < 1, the algorithm tries to exploit the neighboring solutions during the exploitation stages. Essentially, when |P| > 1, exploration occurs while when |P| < 1 exploitation takes place.

3.3 Exploitation phase

The hawks assault the prey target they spotted in the earlier phase during this phase. Prey, however, frequently tries to flee from dangerous situations. According to the prey escape behaviors and pursuit strategies of hawks, four possible hunting strategies emerge. Preys are always attempting to flee dangerous circumstances. Let l be the likelihood of a target successfully fleeing (l < 0.5) or unsuccessfully escaping (l ≥ 0.5) before the attack. Hawks will use a hard siege or a soft siege no matter what the prey does. The parameter P is used to model this method. Thus, soft siege happens for |P| ≥ 0.5 and strong siege occurs for |P| < 0.5 [23].

3.3.4 Hard siege with progressive rapid dives

When l < 0.5 and |P| < 0.5, the hare is exhausted, so the hawks perform a hard siege before the sudden attack. The hare situation is similar to that of the mild siege (Equation (35)), but on this occasion the hawks try to decrease the distance between their position and the hare’s position. While Equation (35) is the same, X and Z now obey new rules:

(36) X = Y prey k − P G ⋅ Y prey k − Y a k

(37) Z = X + U × L evy F S

where Y _a (k) results from Equation (27) [28]. Based on the above, the flowchart of the HHO algorithm is shown in Figure 4, while the corresponding pseudocode is given in Figure 5.

(38) X = Y prey ( k ) − P G ⋅ Y prey ( k ) − Y ( k )

Figure 4:

Flowchart of the HHO algorithm used here.

Figure 5:

Pseudocode of the HHO algorithm used here [28].

4.1 Scenarios examined

In order to examine the system dynamics in terms of transient performance parameters like peak undershoot (U _sh ), peak overshoot (O _sh ), and settling time (T _s ) of the corresponding frequency fluctuations due to a load disturbance, the power system is appropriately modeled and the control strategies are implemented in MATLAB and Simulink environments. The power system model under consideration was used to run a variety of scenarios as shown in Table 1 and the simulation results from each instance are reported next. Each time, the four different bio-inspired algorithms (PSO, GWO, WOA and HHO) were used to evaluate the system’s performance. The same procedure was followed for PI, PID and FOPID controllers. According to Table 1 the following apply:

1. In scenario 1 a dedicated controller is used for each source, the steam unit is considered reheated and the load disturbance in the area (ΔP _L ) is equal to 1 %.

2. In scenario 2, all the above applies except that the load disturbance is equal to 2 %.

For the rest of the scenarios, the load disturbance (ΔP _L ) is kept equal to 1 %.

1. Scenario 3 has the characteristics of Scenario 1 and additionally, GRC as well as GDB are considered for the steam unit.

2. In scenario 4, boiler dynamics are additionally considered for the steam unit.

3. In scenario 5, GRC for the hydro unit is taken into account.

4. Finally, Scenario 6 has the characteristics of Scenario 1, except that a common controller was used for all three sources.

It should be noted that the latter scenario is examined for “negative” demonstration purposes i.e. it will be deducted from the results that common controllers (which may be met in some works) are not advised, whereas dedicated controllers for each source are preferable in LFC studies.

4.2 Results and discussion

Signals are sampled at 50 ms which is considered a very safe sampling frequency as it is 8–10 times larger than most common industrial sampling equipment while the system’s frequency deviation (Δf) time responses are simulated. The simulations have yielded significant results, which are summarized in Figures 7 –9 and Tables 2 –4. With regard to the PI controller, Table 2 shows that the ITAE for the different scenarios and for the different algorithms are (in descending order) the PSO (deterioration on average 64.4 %), WOA (deterioration on average 38.43 %), GWO (deterioration on average 37.39 %) compared to the HHO. Figure 7 provides a visual representation of the frequency deviation response of the PI controller exhibited by the different optimization algorithms. The figure clearly shows that the HHO produced the smallest undershoot in all examined cases, while the PSO produced the highest overshoot. Notably, the HHO exhibited the shortest settling time and, on average, the PSO reaches stability more slowly. It is also observed that as nonlinearities are added to the system (i.e. scenarios 3, 4), the settling time increases in the cases of GWO, WOA and PSO and decreases in the case of HHO.

Figure 7:

System’s frequency deviation time responses using PI controller tuned by PSO, GWO, WOA and HHO, for the different examined operating scenarios (with respect to Table 1): (a) scenario 1, (b) scenario 2, (c) scenario 3, (d) scenario 4, (e) scenario 5 and (f) scenario 6.

Figure 8:

System’s frequency deviation time responses using PID controller tuned by PSO, GWO, WOA and HHO, for the different examined operating scenarios (with respect to Table 1): (a) scenario 1, (b) scenario 2, (c) scenario 3, (d) scenario 4, (e) scenario 5 and (f) scenario 6.

Figure 9:

System’s frequency deviation time responses using FOPID controller tuned by PSO, GWO, WOA and HHO, for the different examined operating scenarios (with respect to Table 1): (a) scenario 1, (b) scenario 2, (c) scenario 3, (d) scenario 4, (e) scenario 5 and (f) scenario 6.

With regard to the PID controller, Table 3 shows that the ITAE for the different scenarios and for the different algorithms are (in descending order) the PSO (deterioration on average 68.2 %), WOA (deterioration on average 47.2 %), GWO (deterioration on average 27.4 %) compared to the HHO. Additionally, when comparing the objective functions of the PI and PID controllers, it is evident that the PID controller consistently achieves a lower ITAE in all examined scenarios. This improvement can be attributed to the presence of the differential term in the PID controller, which significantly reduces the settling time. In addition, by observing the Figures 7 and 8, it can be seen that in the case of a PID controller, the undershoot, overshoot and settling time have been significantly reduced. As shown by the PI controller and similarly observed with the PID controller, the maximum settling time in all scenarios is given by PSO and the minimum by HHO.

For the last controller (FOPID) considered for the LFC problem in the examined system, the ITAE value for the different algorithms (in descending order) was obtained to be the PSO (deterioration on average 80.5 %), WOA (deterioration on average 60.5 %), GWO (deterioration on average 47.8 %) compared to the HHO. In the scenarios with non-linearities HHO-FOPID succeeded in giving lower objective function values than HHO-PID: 30.7 %, 26.7 %, 73.2 %, for Scenarios 3, 4 and 5, respectively. The FOPID controller exhibits an additional distinction not present in the other two controllers. In scenario 6, during the optimization process, the algorithms successfully determined distinct gains for each controller, with HHO yielding the most optimal results.

All algorithms successfully found gain values that led the system to stability for all types of controllers. It was observed that the HHO consistently produced the smallest objective function values across all three controllers. Additionally, indicatively for Scenario 5 — which incorporates all the non-linearities — the frequency deviation response of the HHO-FOPID controller showed a shorter settling time (improved by 75.88 % and 6.49 % for the HHO-PI and HHO-PID controllers, respectively) and less overshoot (improved by 87.72 % and 78.13 % for the HHO-PI and HHO-PID controllers, respectively), as depicted in Figures 7 –9. These findings indicate that the HHO algorithm has the potential to be a valuable tool for optimizing gain values for various controllers. It can enhance system performance and stability. However, it is important to acknowledge the challenge of determining the superiority of one algorithm over another, as the no-free-lunch (NFL) theorem clarifies. This theorem states that no single algorithm can effectively solve every problem. Instead, different algorithms are better suited for different problems. In the context of this study, the HHO algorithm demonstrates superiority in solving the single-area multi-source LFC problem, but it may encounter challenges when applied to other types of problems.

5.1 Random load disturbance

In this examined case, the dynamical load fluctuation of ΔP _L is considered to demonstrate the effectiveness of the proposed HHO-FOPID control method. To show the performance robustness of the HHO-FOPID controller at different perturbations compared with other algorithms, a random load disturbance is applied, as shown in Figure 10. The system response oscillations are illustrated in Figure 11 which depicts the frequency deviation response comparison of the three different scenarios (i.e. scenarios 1, 4 and 5) with the HHO-FOPID controller and indicatively the frequency response for Scenario 1 with the four different bio-inspired algorithms.

Figure 10:

Random load profile.

Figure 11:

Frequency deviation response for system load variation as per Figure 10, (a) with the FOPID controller and for scenario 1, (b) with the HHO-FOPID controller for scenarios 1, 4 and 5.

In this case, the frequency deviation response range lied between ±5 mHz for the PSO algorithm, ±3 mHz for WOA, while the corresponding range for the rest of the algorithms was only ±2 mHz. It can be seen that the proposed HHO-FOPID controller continues to provide the best response for random load variation. It provides less steady-state error, faster response, and better damped oscillation than that achieved using other algorithms. By observing Figure 11(a), it can be seen that, at the instant where the maximum overshoot (t = 80 s) and maximum undershoot (t = 40 s) occurred, the load went from the maximum value to the minimum and from the minimum to the maximum, respectively. The remaining peaks of the frequency deviation response appeared whenever the load changed between the other levels.

Based on the analysis of Figure 11(b), it can be inferred that the frequency deviation response range for scenario 1 is ±1 mHz, whereas for the other two scenarios involving non-linearities, the corresponding range is ±6 mHz. Additionally, it was seen that there were “sharp” edges in this plot in scenarios 4 and 5. This can be attributed to the increased complexity of the systems in these Scenarios due to the presence of the GRC, GDB and boiler. As demonstrated in Figure 11, the HHO-FOPID controller was capable of balancing demand and generation even under fast-changing random load conditions.

5.2 Sensitivity analysis

Sensitivity analysis is crucial for proving any controller’s effectiveness. In this subsection, a sensitivity analysis to assess the stability of the suggested approach (HHO-FOPID) is performed. The effects of the variation of the system constants on the dynamic response of the power system are investigated and the obtained results are summarized in Table 6. Each system constant is increased by 20 % of its nominal value, and the effect on system dynamics is examined in terms of transient performance parameters (i.e. peak undershoot U _sh , peak overshoot O _sh and settling time T _s . It is demonstrated that the impact of modifying system constants on the system’s dynamic performance is almost negligible. It can be concluded that the suggested approach is stable under a variety of scenarios and that the controller gains determined under nominal conditions do not need to be changed when a large range of variation is applied to the system constants.