# Design of HHO-PID Controllers for Load Angle of Power Plant Synchronous Generators

Authors: Morteza Abdolhosseini; Rohollah Abdollahi
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

Accurate analysis of load angle stability of the synchronous generator depends on considering mechanical parameters in addition to electrical parameters in synchronous generator modeling, which leads to increasing state-space variables and model complexity and ultimately reduces the accuracy and performance of conventional controllers. Therefore, an intelligent controller is required for the stable and reliable operation of the synchronous load angle generator. With this approach in this paper, first considering the electrical and mechanical parameters, a complete 14-order model of the synchronous generator is obtained. Then, a PID controller based on the Harris optimization algorithm is proposed to stabilize the load angle of the synchronous generator. The HHO optimization algorithm is a population-based optimization technique that optimizes PID controller coefficients. The effectiveness of the new HHO technique in optimizing the controller PID coefficients over known optimization techniques such as GA-PID and Modified GA-PID has been confirmed using simulation results. The simulation results show that the proposed controller has a higher performance in load angle stability compared to the sliding mode and fuzzy-PID controllers.


---

## Body

## 1. Introduction

Stability is one of the essential characteristics and requirements in dynamic systems. Sustainability has different definitions; these include Lyapunov stability and input-output stability. In physical matters, when we say that the system is unstable, the system’s output goes beyond an acceptable level so that it is not possible to return to the original state without human intervention. For example, an inverting pendulum goes out of balance with the slightest impact and falls to the ground, or an asynchronous generator is unstable when it goes out of sync.If a three-phase short circuit occurs in a power grid and the relays do not work, all the generators will be out of sync mode one after the other. Theoretically, if the physical constraints are removed, the speed of the rotors will increase to infinity. However, due to physical limitations, if the speed is too high, the bearings will be damaged, and the rotor may be dislodged. Although generators and turbines have compelling protection systems that never allow this to happen, the power system is stable when control systems fail to resolve the disturbance, and, to maintain the health of the system, protection systems enter the circuit [1].The main variable in the calculations related to the stability of the generator is the load angle. By increasing the load angle too much, the generator is in danger of exceeding the stable operating range. So, the load angle should be reduced using control methods. But, by reducing the load angle too much, the active output power of the generator will be less than the allowable limit. Therefore, to improve the stable operation of the generator, it is necessary to control the load angle of the generator in such a way that, firstly, the constraints and limits of stability are observed; secondly, the active output power of the generator is optimally adjusted [2]. Power angle or load angle is the angle between the voltage induced by the rotor on the stator and the output voltage of the generator. Another definition of load angle is the difference between the angular position of the rotor and the axis that rotates at synchronous speed. Load angle is the difference in mechanical angle between the magnetic field of the rotor and the magnetic flux of the armature of the synchronous generator [3].In [4], to determine the load angle of the simultaneous prominent pole generator, an optical encoder has been used to detect the position of the rotor and estimate the load angle based on the measured electrical values. Load angle estimation using voltage-current vector diagram and parameters of generators, transformers, and transmission lines in [5] and in [6] based on the generator terminal and measuring the phase angle of the terminal relative to the infinite busbar of the phase measurement unit (PMU) is provided for synchronous generators. Considering the field voltage as input and rotor angle and active power as output, the third-order nonlinear state-space model of the generator is simultaneously identified [7], and the unknown parameters based on the least-squares error method are met [8].In [9], a tool for measuring load angle using real-time digital signal processing (DSP) for synchronous generators is proposed. A comparison of methods for simultaneously defining machine load angles is given in [10]. A load angle controller is installed in the synchronous generator excitation control system. The task of the load controller is to return the generator to a stable operating area, which in the event of a breakdown will lead the generator to stability [11].In [12–14], to adjust the load frequency characteristic, both control structures based on scheduled BFOA-PI and GA-PI, two PI-GSA controllers and adaptive MPC based on bat algorithm, are presented, respectively. To control the voltage regulation, in [15], a constrained genetic algorithm is applied on a nonfragile PID controller, and in [16], a robust PID based on the FSA optimization algorithm is proposed. In [17], to control the blade pitch in wind energy systems, an adaptive MPC controller based on the crow search algorithm is used. A hybrid of fuzzy logic and nonlinear sliding mode control to control the distributed controllable loads is suggested in [18].In [19], modification of control objective is proposed to reduce the rotor angle acceleration power systems and improve the stability margin during system disturbances. In [20], the design of a proportional derivative controller (PID) based on low-frequency oscillation damping using a power system stabilizer for large-scale photovoltaic farms is proposed. Also, a method based on optimal controller coordination based fractional-order proportional-integral-derivative for low-frequency oscillation damping in the power system is presented in [21]. In [22], a fractional-order nonlinear proportional-integral-derivative (FONPID) controller for wind energy system is discussed. A proportional-integral-derivative based on reinforced learning neural network for discrete mode load frequency control problems in hydrothermal hybrid distributed generation has been studied in [23]. In [24], to improve the voltage and reactive power, a particle swarm optimization based on a hybrid bacterial foraging optimization algorithm is proposed to optimize the control of STATCOM and AVR.Increasing the load on a generator will increase the load angle, and raising more than allowed will bring the generator operating point into the unstable working area. In many power systems, no device is responsible for controlling the load angle independently. However, at operating points that are far enough away from the stability range, load angle fluctuations are damped. The problem arises when a large enough disturbance is applied to the system, and the system enters the unstable zone and loses synchronicity by increasing the load angle too much. Under such circumstances, having a controller that reduces the load angle to the allowable value can prevent the system from becoming unstable and out of sync. Therefore, researchers have always been looking for solutions to control the load angle of the generator and safer operation of the power system.Table1 shows the classification of the load angle identification and control methods. As can be seen, the proposed control method uses the power plant synchronous generator with 14th order model, which is higher in rank than the existing methods and therefore more complex and accurate. The simulation results confirm the high performance of the proposed control system in this complete model.Table 1
Classification of the load angle identification and control methods.
Ref.Load angle estimationGenerator modelingLoad angle controlHardware implementation and real-time checkOptical encoderVoltage-current vector diagramPhasor measurement units (PMUs)Least squares errorSys. Order-3Sys. Order-7Sys. Order-14[4]●[5]●[6]●[7]●[8]●[9]●[10]●[11]●●[28, 29]●●Proposed method●●In [25], a reduced-order aggregate model based on the balanced truncation is proposed in DC microgrid. Also, in [26], a closed-loop transfer function model with reduced order is proposed for the real-time simulation of power systems. However, the models of modern power plant generator systems are more and more complex, the operating conditions are complex and variable, and different nonlinear factors are paired with each other. As a result, many model reduction methods are no longer applicable, and, by reducing the order of the model, the effect of all mechanical variables or all electrical variables may be ignored. This is a challenge to study the dynamic properties of complex high-dimensional systems and to design optimized system parameters. In addition, special attention should be paid to reducing the model of interconnected systems in order to maintain the integrity and connection structure between different subsystems. Therefore, using the HHO-PID method in this type of system will be challenging and may have an error in the controller output and not be accurate enough [27].With this approach, in this paper, a comprehensive and complex model with 14 variables has been selected and the HHO-PID control method has been used to control the load angle by considering 8 electrical and 6 mechanical variables simultaneously. The major contributions of the paper include the following:(i)
A complete 14-order model of the synchronous generator, with simultaneous consideration of electrical and mechanical parameters, has been realized to accurately investigate the load angle of the synchronous generator(ii)
A PID controller based on the HHO algorithm is proposed to ensure the stability of the load angle of the synchronous generator in the 14-order model(iii)
The HHO optimization algorithm is a population-based optimization technique that optimizes PID controller coefficients(iv)
The effectiveness of the proposed algorithm compared to conventional intelligent algorithms such as GA and Modified GA as well as the sliding mode and fuzzy-PID controllers has been verified using software analysisThe rest of this article is arranged as follows: In Section2, the equations of the state space of a simple power plant are expressed, and, in Sections 3 and 4, the PID controller and HHO algorithm are introduced. In Section 5, the simulation results are presented. Section 6 concludes the study.


## 2. State-Space Equations of a Simple Power Plant

This section examines a complete power plant model used in many power system dynamics studies. This system is connected to an infinite busbar through a transformer and transmission lines.The differential equations for this model are(1)δ˙=Δω,Δω˙=w0Tm−Te−KdΔω2H,ψ˙fd=w0Vfd−RfdIfd,ψ˙d=w0Vbd+ψq−IdRa+Re+ψqΔω,ψ˙kd=−w0RkdIkd,ψ˙q=w0Vbq−ψd−IqRa+Re−ψdΔω,ψ˙kq=−w0RkqIkq.The excitation system model of this generator is a first-order system, which represents a static exciter. Its distinctive feature is a small-time constant and limiting the exciter output to both positive and negative limits. The static excitation system equation is also considered as a first-order equation as follows:(2)V˙fd=VR−VfdSex.Total mechanical torque (Tm) is obtained by combining the torques produced in high-, medium-, and low-pressure turbines, related to the FHP, FIP, and FLP coefficients in the model in Figure 1. When the variables are calculated per unit, the sum of these coefficients is one. Figure 1 shows a nonlinear model of a turbine suitable for transient dynamic studies.Figure 1
Nonlinear model of the turbine.In a mechanical power generating system, the necessary steam is first generated by the boiler. It then enters the high-pressure (H.P) class of the turbine with the main valve. The steam from the high-pressure floor enters the reheater and then the medium-pressure (I.P) floor of the turbine, where the amount of steam transferred between the two floors is also controlled by the shut-off valve. Finally, the power turbine will provide three power levels of high pressure, medium pressure, and low pressure. The differential equations for the different parts of the turbine and the reheater are also presented as follows:(3)Y˙HP=GVMP0−YHPSHP,Y˙RH=YHP−YRHSRH,Y˙IP=GVIYRH−YIPSIP,Y˙LP=YIP−YLPSLP,G˙VM=UGM−GVMSGVM,G˙VI=UGI−GVISGVI.It should be noted that the total impedance includes the impedance of the transformer and the transmission line connected to it in constant values according to the following equations:(4)Re=Rtr+RL,xe=xtr+xL.The infinite bus voltage is also a constant value per unit.The following equations also apply to the electrical torqueTe and the mechanical torque Tm.(5)Te=ψdId−ψqId,Tm=FHPYHP+FIPYIP+FLPYLP.Figure2 shows the block diagram of electric torque generation. In this block diagram, according to equation (1), ψq and ψd are obtained first and then they are used to generate electric torque.Figure 2
Block diagram of electric torque generation.Now, by examining the details of this model, extracting some other relations, and combining them in the form of first-order differential equations, the nonlinear state-space equations are obtained for the system under study. The first step is to determine the state variables or, in other words, to determine the state vector of the system and the system input vector. In this 14-order synchronous generator, two orders are associated with mechanical modes and rotor movement, five other orders are related to the flow of stator, damper, and exciter windings, one order is related to generator excitation system, and six degrees are dependent on mechanical power generating system.Thus, we have the following vector:(6)X=δ,Δω,ψfd,ψd,ψkd,ψq,ψkq,Vfd,YHP,YRH,YIP,YLP,GVM,GVI.It represents the vector of system state-space variables with dimensions of 14 × 1 and there is the following vector:(7)U=UgVRT,Ug=UGM=UGI.It also represents the system input vector with dimensions of 2 × 1.The next step is to determine expressionsId, Iq, Ikd, Ikd, and Ifd in terms of the state-space variables defined in equation (6). The differential equations including these state variables can be generalized to the state-space equations. Now, according to the definition of the state-space variables and the desired system inputs after simplification and sorting, the nonlinear state-space equations for this synchronous generator can be obtained, which are as follows:(8)x˙1=x2,x˙2=w0FHPx9+FIPx11+FLPx12−x6y1dx4+y4dx3+y5dx5+x4y1qx6+y3qx7−Kdx22H,x˙3=w0x8−Rfdy4dx4+y2dx3+y6dx5,x˙4=w0Vbsinx1+x6−Ra+Rey1dx4+y4dx3+y5dx5+x2x6,x˙5=−w0Rkdy5dx4+y6dx3+y3dx5,x˙6=w0Vbcosx1−x4−Ra+Rey1qx6+y3qx7−x2x3,x˙7=−w0Rkqy3qx6+y2qx7,x˙8=u2−x8Sex,x˙9=P0x13−x9SHP,x˙10=x9−x10SRH,x˙11=x10x14−x11SIP,x˙12=x11−x12SLP,x˙13=u1−x13SGVM,x˙14=u1−x14SGVI.Equation (8) can be represented by the closed form of nonlinear equations of equation (9). In equation (9), the rows of matrix f are the same as the functions of equation (8).(9)x˙=fx.u.To linearize this model around a given operating pointx0,u0 in the next n+m space, the dimensions of matrices A and B will be n×n and n×m, respectively. Now, to convert the linear state-space model, the condition fx0,u0=0 must be met. In other words, vector x0 is the same as the permanent state vector of the system, which indicates that the system has reached its steady state and the stagnation condition x˙=0 is met at the equilibrium point. So, to determine each x0,u0 equilibrium point, the optimal approach is as follows:(10)fx0,u0=0.To linearize a synchronous machine model, we must first calculate the load angle,id,iq,vd, and vq, according to the following equation:(11)Rt≜Ra+Re,xt≜xq+xe,I≜Pt2+Qt2Vb,φ≜tan−1QtPt,x10=tan−1xtIcosφ−RtIsinφVb+RtIcosφ+xtIsinφ,vbd=Vbsinx10,vbq=Vbcosx10,Id=Isinx10+φ,Iq=Icosx10+φ.In equation (11), Pt, Qt, and Vb represent the preunit values of active power, reactive power, and infinite bus voltage, respectively, which are used as known factors in determining the operating point of the system’s steady state. Therefore, for specific values of these factors, 7 permanent state values for the state variables x2 to x8, which are represented by the symbols x20 to x80, respectively, are calculated as follows:(12)x20=0,x40=Vbq+Iq×Rt,Ifd=x40+IdXd+XeXad,x30=Xfd×Ifd−Xad×IdIfd,x50=Xad×Ifd−Xad×Id,x60=−Vbd+Id×Rt,x70=−Xaq×Iq,x80=Rfd×Ifd.Since there is a multiplication of two state variables in the governor model, a quadratic equation is obtained to calculate the governor’s input value in the steady state:(13)P0FLP+P0FIPu120+P0FHPu10+x40y1qx60+x40y3qx70−x30y4dx60−x50y5dy60−x40y1dx60=0,which is a quadratic algebraic equation in terms of the unknown variable u10. It can be used to determine the value of the steady-state input u1 at the equilibrium point of the system (one of the roots is acceptable). In addition, other equilibrium points are obtained from the following equation:(14)u20=x80,x90=x100=P0u10,x110=x120=P0u120,x130=x140=u10.Finally, the linearized model will bex˙n×1=Ax+Bum×1. As shown in the above equations, for each given pair Pt,Qt with the initial default Vb=1p.u the coordinates of the equilibrium vector x0 and the input vector u0 are determined.The steady-state values of the state variables are given as follows:(15)xi0=0.7802,0,1.5562,0.7696,1.3645,−0.8071,−0.6053,0.0026,1.0430,1.0430,1.0879,,1.0879,1.0430,1.0430.By linearizing around the operating point, the system transfer function is obtained as follows:(16)system=19.41s10+2372s9+4.513e04s8+2.07e05s7+3.301e05s6+2.666e05s5+2.882e05s4+6.198e04s3+3177s2+37.2s+0.06977s14+218.2s13+1.375e04s12+2.036e05s11+1.103e06s10+2.624e06s9+3.106e06s8+2.889e06s7+2.178e06s6+4.897e05s5+1.62e05s4+2.144e04s3+923.4s2+10.29s+0.01911.


## 3. Introducing HHO Optimization Algorithm

In this paper, a population-based and nature-inspired Harris Hawk Optimization (HHO) algorithm is used. The main idea of the HHO algorithm is the participatory behavior and pursuit style of Harris hawk in nature, which is known as a surprise attack [30].In this intelligent strategy, several suspicious hawks are tried to surprise prey from different directions.The Harris hawk can exhibit various chase patterns based on the dynamic nature of hunting scenarios and routines. For this purpose, a mathematical model is proposed, and then a stochastic metaheuristic algorithm based on the proposed mathematical model is designed to deal with various optimization problems [31].In this section, we model the exploration and operation phases of the HHO algorithm by exploring prey, surprise pounce, and different attacking strategies of Harris hawks.The HHO algorithm is a gradient-independent, population-based optimization technique, so it can be applied to solve any optimization problem provided there is a suitable formula.


### 3.1. Exploration Phase

In this section, we explain the mechanism of HHO exploration. If we consider the nature of the Harris hawk, it can identify and track prey using strong eyes, but occasionally prey cannot be easily seen. As a result, the hawk waits, observing the desert to identify prey after several hours.In the HHO algorithm, hawks are considered candidate solutions, and the best solution at each stage is considered the desired prey or almost optimal answer. In HHO, hawks are randomly placed in different areas, waiting to identify targets based on two strategies. We consider the same chance for each of the strategies.Strategy 1.
They determine their position based on the position of other family members and the position of the prey (rabbit) (to be close enough to them at the time of the attack), which is formulated in equation (18) for the condition q<0.5.Strategy 2.
Hawks are randomly placed on tall trees (random locations around the group), which is modeled in equation (18) with condition q≥0.5.(17)Xt+1=Xrandt−r1Xrandt−2r2Xtq≥0.5,Xrabbitt−Xmt−r3LB+r4UB−LB,q<0.5.
In equation (17), Xt+1 is the position of the hawks in the next iteration t, Xrabbitt is the position of the rabbit, Xt is the current position of the hawks, and the coefficients r1, r2, r3, r4, and q are random numbers in the range 0,1, which are updated in each iteration. The LB and UB parameters show the lower and upper limits of the variables, respectively. Xrandt is a hawk randomly selected from the current population, and Xmt is the average position of the current population of hawks.
A simple model for random generation of locations within the range (LB, UB) is considered. In the first strategy, solutions are generated based on random location and other falcons. In the second strategy, the difference between the best location so far and the average group position plus a random scale component based on the range of variables is considered.
The average position of the hawks is calculated with the following equation:(18)Xmt=1N∑i=1NXit,where Xit indicates the position of each hawk in the iteration of t and N denotes the total number of hawks.


### 3.2. Transition from the Exploration Phase to the Exploitation Phase

The HHO algorithm can be switched from the exploration phase to exploitation and then switch between different exploitation behaviors based on the residual energy of the prey. The energy of the prey decreases over time as it escapes. Prey energy is modeled based on the following equation:(19)E=2E01−tT,where E is the escape energy of the prey, T is the maximum number of iterations of the algorithm, and E0 is the initial amount of energy.In HHO, the value ofE0 in each iteration is randomly selected in the range from −1 to 1.When the value ofE0 decreases from 0 to −1, the rabbit becomes physically weak, and when the value of E0 increases from 0 to 1, it means that the rabbit becomes stronger. Dynamic escape energy E has a decreasing trend during the iterations of the algorithm. When the escape energy is E≥1, the falcons search different areas to locate the rabbit, so HHO will be in the exploration phase and when E<1 it will be in the exploitation phase.


### 3.3. Exploitation Phase

At this stage, the hawks perform the surprise attack behavior by attacking the prey detected in the previous step. But they usually try to escape from dangerous situations. As a result, in real cases, different chasing behaviors are formed.Based on prey escape behaviors and hawks chasing strategies, four possible approaches have been proposed in HHO to model the exploitation process.Prey always tries to escape from threatening situations. Assume that ther prey has a chance of a successful escape (r<0.5) or a failed escape (r≥0.5) before the surprise pounce.As the prey runs, the hawks form a hard or soft siege based on their residual energy to catch the hunt. In real life, the hawks get closer and closer to the prey to increase their group chance of killing the rabbit by surprise pounce. After a few minutes, the rabbit gradually loses its energy, and then the hawks intensify the siege process to catch the tired prey easily.To model the above strategy, parameterE is defined. In this case, a soft surround occurs when E≥0.5, and a hard surround occurs when E<0.5.


#### 3.3.1. Soft Siege

WhenE≥0.5 and r≥0.5, the rabbit still has enough energy and tries to escape with misleading jumps but ultimately cannot. During these efforts, the hawks gently surround him to make the rabbit more tired and then launch a surprise pounce.This behavior is modeled according to the following equation:(20)Xt+1=ΔXt−EJXrabbitt−Xt,ΔXt=Xrabbitt−Xt,where ΔXt is the difference between the rabbit position vector and the current position in the iteration of t, r5 is a random number in the range of 0-1, and J=21−r5 indicates the strength of the rabbit’s unexpected jump during the escape process. The value of J in each iteration changes randomly to simulate the nature of rabbit movement.


#### 3.3.2. Hard Siege

Whenr≥0.5 and E≤1, the prey is tired enough, and its escape energy is reduced. In addition, the hawks severely encircle the prey to make a surprise pounce. In this case, the current position is updated based on the following equation:(21)Xt+1=Xrabbitt−EΔXt


#### 3.3.3. Soft Siege with Fast Forward Attacks

WhenE≥0.5 but r<0.5, the rabbit has enough energy to escape, and the soft siege will still be in place. To mathematically model prey escape patterns and mutation motions, the concept of Levy or LF flight is used in the HHO algorithm.The LF concept is used for the deceptive movements of prey zigzags (especially rabbits) during the escape phase, as well as the irregular and rapid movements of hawks around the escaped prey. Hawks make several team immediate drives around the rabbit and gradually correct their location and path due to the prey’s deceptive movements. Inspired by the actual behavior of falcons, it is thought that they can choose the best position progressively to attack prey when they want to catch prey in a competitive environment.Therefore, to perform a soft siege, we assumed that the hawks could evaluate (decide) their next move based on the following equation:(22)Y=Xrabbitt−EJXrabbitt−Xt.Then, they compare the probable outcome of such a move with the previous dive to see if it would be a good dive. When they see the prey making more deceptive movements, when they approach the rabbit, they also start to make irregular, sudden, and fast dives. They are thought to dive into LF-based patterns using the following equation:(23)Z=Y+S×LFD,where D is the dimension of the problem, S is a random vector of size D, and LF is a Levy flight function, calculated by the following equation:(24)LFx=0.01×u×σv1/β,σ=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β,where u and v are random numbers in the range of 0-1 and β is a constant value equal to 1.5.Hence, the final strategy for updating the positions of the hawks in the soft siege stage is done by the following equation:(25)Xt+1=YifFY<FXt,ZifFZ<FXt,where Y and Z are determined by equations (22) and (23).In each step, only the best position betweenY or Z is selected as the next location. This strategy applies to all search agents.


#### 3.3.4. Hard Siege with Fast Forward Attacks

When |E| < 0.5 and r < 0.5, rabbits do not have enough energy to escape, and a hard siege is carried out before a surprise pounce to catch and kill prey. The position of this stage on the prey side is similar to the situation in the soft siege, but the hawks try to reduce their average distance with the fleeing prey. Therefore, equation (26) is performed in the event of a hard siege:(26)Xt+1=YifFY<FXt,ZifFZ<FXt,where Y and Z are obtained as in the following equation:(27)Y=Xrabbitt−EJXrabbitt−Xmt,Z=Y+S×LFD.We have that Xmt is calculated using equation (18).


## 3.1. Exploration Phase

In this section, we explain the mechanism of HHO exploration. If we consider the nature of the Harris hawk, it can identify and track prey using strong eyes, but occasionally prey cannot be easily seen. As a result, the hawk waits, observing the desert to identify prey after several hours.In the HHO algorithm, hawks are considered candidate solutions, and the best solution at each stage is considered the desired prey or almost optimal answer. In HHO, hawks are randomly placed in different areas, waiting to identify targets based on two strategies. We consider the same chance for each of the strategies.Strategy 1.
They determine their position based on the position of other family members and the position of the prey (rabbit) (to be close enough to them at the time of the attack), which is formulated in equation (18) for the condition q<0.5.Strategy 2.
Hawks are randomly placed on tall trees (random locations around the group), which is modeled in equation (18) with condition q≥0.5.(17)Xt+1=Xrandt−r1Xrandt−2r2Xtq≥0.5,Xrabbitt−Xmt−r3LB+r4UB−LB,q<0.5.
In equation (17), Xt+1 is the position of the hawks in the next iteration t, Xrabbitt is the position of the rabbit, Xt is the current position of the hawks, and the coefficients r1, r2, r3, r4, and q are random numbers in the range 0,1, which are updated in each iteration. The LB and UB parameters show the lower and upper limits of the variables, respectively. Xrandt is a hawk randomly selected from the current population, and Xmt is the average position of the current population of hawks.
A simple model for random generation of locations within the range (LB, UB) is considered. In the first strategy, solutions are generated based on random location and other falcons. In the second strategy, the difference between the best location so far and the average group position plus a random scale component based on the range of variables is considered.
The average position of the hawks is calculated with the following equation:(18)Xmt=1N∑i=1NXit,where Xit indicates the position of each hawk in the iteration of t and N denotes the total number of hawks.


## 3.2. Transition from the Exploration Phase to the Exploitation Phase

The HHO algorithm can be switched from the exploration phase to exploitation and then switch between different exploitation behaviors based on the residual energy of the prey. The energy of the prey decreases over time as it escapes. Prey energy is modeled based on the following equation:(19)E=2E01−tT,where E is the escape energy of the prey, T is the maximum number of iterations of the algorithm, and E0 is the initial amount of energy.In HHO, the value ofE0 in each iteration is randomly selected in the range from −1 to 1.When the value ofE0 decreases from 0 to −1, the rabbit becomes physically weak, and when the value of E0 increases from 0 to 1, it means that the rabbit becomes stronger. Dynamic escape energy E has a decreasing trend during the iterations of the algorithm. When the escape energy is E≥1, the falcons search different areas to locate the rabbit, so HHO will be in the exploration phase and when E<1 it will be in the exploitation phase.


## 3.3. Exploitation Phase

At this stage, the hawks perform the surprise attack behavior by attacking the prey detected in the previous step. But they usually try to escape from dangerous situations. As a result, in real cases, different chasing behaviors are formed.Based on prey escape behaviors and hawks chasing strategies, four possible approaches have been proposed in HHO to model the exploitation process.Prey always tries to escape from threatening situations. Assume that ther prey has a chance of a successful escape (r<0.5) or a failed escape (r≥0.5) before the surprise pounce.As the prey runs, the hawks form a hard or soft siege based on their residual energy to catch the hunt. In real life, the hawks get closer and closer to the prey to increase their group chance of killing the rabbit by surprise pounce. After a few minutes, the rabbit gradually loses its energy, and then the hawks intensify the siege process to catch the tired prey easily.To model the above strategy, parameterE is defined. In this case, a soft surround occurs when E≥0.5, and a hard surround occurs when E<0.5.


### 3.3.1. Soft Siege

WhenE≥0.5 and r≥0.5, the rabbit still has enough energy and tries to escape with misleading jumps but ultimately cannot. During these efforts, the hawks gently surround him to make the rabbit more tired and then launch a surprise pounce.This behavior is modeled according to the following equation:(20)Xt+1=ΔXt−EJXrabbitt−Xt,ΔXt=Xrabbitt−Xt,where ΔXt is the difference between the rabbit position vector and the current position in the iteration of t, r5 is a random number in the range of 0-1, and J=21−r5 indicates the strength of the rabbit’s unexpected jump during the escape process. The value of J in each iteration changes randomly to simulate the nature of rabbit movement.


### 3.3.2. Hard Siege

Whenr≥0.5 and E≤1, the prey is tired enough, and its escape energy is reduced. In addition, the hawks severely encircle the prey to make a surprise pounce. In this case, the current position is updated based on the following equation:(21)Xt+1=Xrabbitt−EΔXt


### 3.3.3. Soft Siege with Fast Forward Attacks

WhenE≥0.5 but r<0.5, the rabbit has enough energy to escape, and the soft siege will still be in place. To mathematically model prey escape patterns and mutation motions, the concept of Levy or LF flight is used in the HHO algorithm.The LF concept is used for the deceptive movements of prey zigzags (especially rabbits) during the escape phase, as well as the irregular and rapid movements of hawks around the escaped prey. Hawks make several team immediate drives around the rabbit and gradually correct their location and path due to the prey’s deceptive movements. Inspired by the actual behavior of falcons, it is thought that they can choose the best position progressively to attack prey when they want to catch prey in a competitive environment.Therefore, to perform a soft siege, we assumed that the hawks could evaluate (decide) their next move based on the following equation:(22)Y=Xrabbitt−EJXrabbitt−Xt.Then, they compare the probable outcome of such a move with the previous dive to see if it would be a good dive. When they see the prey making more deceptive movements, when they approach the rabbit, they also start to make irregular, sudden, and fast dives. They are thought to dive into LF-based patterns using the following equation:(23)Z=Y+S×LFD,where D is the dimension of the problem, S is a random vector of size D, and LF is a Levy flight function, calculated by the following equation:(24)LFx=0.01×u×σv1/β,σ=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β,where u and v are random numbers in the range of 0-1 and β is a constant value equal to 1.5.Hence, the final strategy for updating the positions of the hawks in the soft siege stage is done by the following equation:(25)Xt+1=YifFY<FXt,ZifFZ<FXt,where Y and Z are determined by equations (22) and (23).In each step, only the best position betweenY or Z is selected as the next location. This strategy applies to all search agents.


### 3.3.4. Hard Siege with Fast Forward Attacks

When |E| < 0.5 and r < 0.5, rabbits do not have enough energy to escape, and a hard siege is carried out before a surprise pounce to catch and kill prey. The position of this stage on the prey side is similar to the situation in the soft siege, but the hawks try to reduce their average distance with the fleeing prey. Therefore, equation (26) is performed in the event of a hard siege:(26)Xt+1=YifFY<FXt,ZifFZ<FXt,where Y and Z are obtained as in the following equation:(27)Y=Xrabbitt−EJXrabbitt−Xmt,Z=Y+S×LFD.We have that Xmt is calculated using equation (18).


## 3.3.1. Soft Siege

WhenE≥0.5 and r≥0.5, the rabbit still has enough energy and tries to escape with misleading jumps but ultimately cannot. During these efforts, the hawks gently surround him to make the rabbit more tired and then launch a surprise pounce.This behavior is modeled according to the following equation:(20)Xt+1=ΔXt−EJXrabbitt−Xt,ΔXt=Xrabbitt−Xt,where ΔXt is the difference between the rabbit position vector and the current position in the iteration of t, r5 is a random number in the range of 0-1, and J=21−r5 indicates the strength of the rabbit’s unexpected jump during the escape process. The value of J in each iteration changes randomly to simulate the nature of rabbit movement.


## 3.3.2. Hard Siege

Whenr≥0.5 and E≤1, the prey is tired enough, and its escape energy is reduced. In addition, the hawks severely encircle the prey to make a surprise pounce. In this case, the current position is updated based on the following equation:(21)Xt+1=Xrabbitt−EΔXt


## 3.3.3. Soft Siege with Fast Forward Attacks

WhenE≥0.5 but r<0.5, the rabbit has enough energy to escape, and the soft siege will still be in place. To mathematically model prey escape patterns and mutation motions, the concept of Levy or LF flight is used in the HHO algorithm.The LF concept is used for the deceptive movements of prey zigzags (especially rabbits) during the escape phase, as well as the irregular and rapid movements of hawks around the escaped prey. Hawks make several team immediate drives around the rabbit and gradually correct their location and path due to the prey’s deceptive movements. Inspired by the actual behavior of falcons, it is thought that they can choose the best position progressively to attack prey when they want to catch prey in a competitive environment.Therefore, to perform a soft siege, we assumed that the hawks could evaluate (decide) their next move based on the following equation:(22)Y=Xrabbitt−EJXrabbitt−Xt.Then, they compare the probable outcome of such a move with the previous dive to see if it would be a good dive. When they see the prey making more deceptive movements, when they approach the rabbit, they also start to make irregular, sudden, and fast dives. They are thought to dive into LF-based patterns using the following equation:(23)Z=Y+S×LFD,where D is the dimension of the problem, S is a random vector of size D, and LF is a Levy flight function, calculated by the following equation:(24)LFx=0.01×u×σv1/β,σ=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β,where u and v are random numbers in the range of 0-1 and β is a constant value equal to 1.5.Hence, the final strategy for updating the positions of the hawks in the soft siege stage is done by the following equation:(25)Xt+1=YifFY<FXt,ZifFZ<FXt,where Y and Z are determined by equations (22) and (23).In each step, only the best position betweenY or Z is selected as the next location. This strategy applies to all search agents.


## 3.3.4. Hard Siege with Fast Forward Attacks

When |E| < 0.5 and r < 0.5, rabbits do not have enough energy to escape, and a hard siege is carried out before a surprise pounce to catch and kill prey. The position of this stage on the prey side is similar to the situation in the soft siege, but the hawks try to reduce their average distance with the fleeing prey. Therefore, equation (26) is performed in the event of a hard siege:(26)Xt+1=YifFY<FXt,ZifFZ<FXt,where Y and Z are obtained as in the following equation:(27)Y=Xrabbitt−EJXrabbitt−Xmt,Z=Y+S×LFD.We have that Xmt is calculated using equation (18).


## 4. Methodology

The first step in defining a model is to introduce its factors. In the first part of this program, the characteristics of the fourteen-degree system model are described. In the second part of the program, the operating points including real power, reactive power, and infinite bus voltage are determined. These values as per unit, 1, 0.5, and 1.05, are considered, respectively. Once the operating point of the system is specified, the values of the steady state can be determined. Then, using them, the values of linear model factors are determined.Now, using the linear model with load angle output, a suitable controller should be designed to improve the load angle characteristics. Due to the performance of the PID controller in the industry, this controller is also used here. The PID controller coefficients are optimized and adjusted by the HHO optimization algorithm. Figure3 shows the flowchart used in this study.Figure 3
Simulation flowchart.The structure of the control ring used is as shown in Figure4. According to this figure, the PID controller is placed before the 14-order generator model, and KP, KI, and KD coefficients of the controller are adjusted and optimized by the HHO algorithm. The objective function that the HHO algorithm minimizes is as follows:(28)Z=w1×OS+w2×ST+w3×RT+w4×SI.Figure 4
Structure of the control loop.In equation (28),w1 to w4 are the coefficients of the objective function, its overshoot OS, ST settling time, and RT rise time, and stability index (SI) is defined as follows:(29)SI=−1minmaxrealpoleT,0,where T is the closed-loop transfer function. For the system to be stable, the real part of all roots of the closed-loop transfer function must be less than zero. In equation (29), first, the largest real value of the roots is calculated and compared to zero. If the root is greater than zero, the SI becomes infinite. Otherwise, the index can have a value close to zero.


## 5. Simulation Results

The parameters of the power plant synchronous generator for simulation are listed in Table2. It should be noted that one generator has been investigated in the simulation. The parameters of the HHO algorithm include the upper and lower limits of the variables, the number of hawks, and the number of iterations as shown in Table 3. The parameters of the HHO algorithm have been taken as a trial and error to obtain better optimization results.Table 2
Simulation parameters of the power plant synchronous generator.
GeneratorTurbine and governorH3.25 puP01.0Ra0.005 puFHP0.24Kd0.025SHP0.3 secRfd0.0015 puSRH10 secxfd1.97 puFIP0.34xad1.86 puSIP0.3 secxaq1.77 puSGVM0.1 secRkd0.0078 puFLP0.42xkd1.94 puSLP0.72 secRkq0.0084 puSGVI0.01 secxkq1.96 puExcitation systemxd2.0 puSex0.01 secxq1.91 puTransformer and transmission linew0=2πf0100π rad/secRtr0.038 puxtr0.1 puRL0.025 puxl0.35 puTable 3
Parameters of HHO algorithm.
Parameters of HHO algorithmNumber of search agentsMaximum number of iterationsNumber of decision variablesLower boundUpper bound100503[1.01, 0.01, 1.01][10, 1, 30]Now in GA-PID algorithm, 1 < KP < 10, 0 < KI < 1, and 1 < KD < 30 are considered, as well as the weights of the objective functionw1 = 0.6, w2 = 0.2, w3 = 0.1, and w4 = 0.1. The reason for choosing the range for the PID controller coefficients is that the best possible solution is obtained for the controller and the optimization algorithm is not trapped in the local minimums. Also, by applying the HHO algorithm, the optimal coefficients KP, KI, and KD of the PID controller are equal to 8, 0.53, and 29.909, respectively.The step response diagrams of the 14 state-space variables affected by the first and second inputs are shown in Figure5. In this figure, Xi are state-space variables and X1 = δ, X2 = Δω, X3 = ψfd, X4 = ψd, X5 = ψkd, X6 = ψq, X7 = ψkq, X8 = Vfd, X9 = YHP, X10 = YRH, X11 = YIP, X12 = YLP, X13 = GVM, and X14 = GVI.Figure 5
Step response diagram of 14 state-space variables.
(a)(b)(c)(d)(e)As can be seen, the step response of the variable statex1 (δ) is affected only by the first input (Ug) and is always zero compared to the change of the second input (VR).As shown in Figure5, only the third, fifth, and eighth state variables have nonzero values for the second input, and the other state variables have zero values for this input. Also, the behavior of the variables of the first, second, fourth, and sixth modes concerning the first input is almost the same.


### 5.1. Active Power Stability with the Proposed Controller

The relationship between output power and load angle is expressed in the following equation [32]:(30)P=Pmaxsinδ.Figures6 and 7 show the output power without the controller and with the controller, respectively. In Figure 6, without applying the controller, the output power has an oscillating and unstable response. With the application of HHO-PID, the output power fluctuations are reduced to zero in less than 10 seconds.Figure 6
Output power before applying the controller.Figure 7
Output power after applying the controller.Figure 8
Nyquist diagram of the load angle.To evaluate the stability of the proposed controller, the Nyquist and Bode diagrams of the load angle are plotted in Figures8 and 9. As shown in Figure 8, this diagram is not circumvented, −1. According to the Nyquist stability criterion, since the open-loop transmission function does not have an instability pole, the load angle output is stable. Using Figure 9, since the gain margin and phase margin are positive, the HHO-PID controller will stabilize the load angle. Figures 8 and 9 clearly show the proposed controller performance in the load angle stability of the generator.Figure 9
Bode diagram of the load angle.Figure 10
Convergence curve of the proposed HHO-PID algorithm.


### 5.2. Proposed Controller Performance

As can be seen in Figure10, after 100 repetitions of the cost function, it has reached 27.1, which confirms the convergence of the proposed method. To evaluate the performance of the proposed controller, the step responses of the HHO-PID controller with fuzzy-PID (FPID) and sliding mode controllers (SMC) have been compared. Therefore, by applying a unit step as Ug input to the synchronous generator system, its load angle characteristic is obtained according to Figure 11(a). As shown in this figure, the characteristic of the load angle is wildly oscillating. Its overshoot and settling time are high, and, after 3000 seconds, it still oscillates, which is not desirable. As shown in Figure 11(e), the step response has a much better performance in its overshoot, settling time, fluctuations, and rise time.Figure 11
Load angle characteristic for a step input. (a) Without the controller, (b) with SMC in control level mode greater than one (SMCU), (c) with fuzzy-PID controller, (d) with SMC in control level mode less than one (SMCD), and (e) with HHO-PID controller.
(a)(b)(c)(d)(e)Also, the performances of the sliding mode controller (SMC) and the fuzzy-PID controller (FPID) are examined. By selecting the control level less than 1 in the Simulink model, the load angle output is obtained as in Figure11(d), which is not a suitable output. Figure 11(b) shows the load angle of a synchronous generator rotor with an SMC in a situation where the control level is greater than one. In this case, it can be seen that this controller does not work properly. Mamdani fuzzy controller with two inputs and three outputs as well as 49 rules is programmed. As shown in Figure 11(c), the load angle is completely unstable and therefore indicates that the FPID controller cannot be an excellent option to control the load angle of the synchronous generator rotor.In SMC and FPID controller, due to the fact that the power plant generator model uses mechanical parameters and the system order is 14, these controllers have not performed well and have diverged. In [28], the load angle of the 7th-order synchronous generator system is controlled by a PID controller using a genetic algorithm. A comparison of the performances of these GA-PID [28], Modified GA-PID, SMCU, SMCD, and FPID methods is presented in Table 4. The results show that the proposed algorithm and system are reduced by 89% at rise time, 72.6% at settling time, and 1.59% at overshoot and have better performance than GA-PID in all step response parameters. In the following, the optimization algorithm used is compared with the modified genetic algorithm. As can be seen in Table 4, the rise time and settling time of the proposed algorithm have been reduced by 36.6% and 74.2% compared to the genetic algorithm, respectively.Table 4
Comparison of the proposed controller performance with GA-PID [28], Modified GA-PID, SMCU, SMCD, and FPID.
Rise time (sec)Settling time (sec)Overshoot (%)Peak time (sec)HHO-PID0.436913.700854.12641.2101GA-PID [28]450558Modified GA-PID0.689753.187852.66781.9944SMCU1500——2000SMCD1022004820FPID2400——2440


## 5.1. Active Power Stability with the Proposed Controller

The relationship between output power and load angle is expressed in the following equation [32]:(30)P=Pmaxsinδ.Figures6 and 7 show the output power without the controller and with the controller, respectively. In Figure 6, without applying the controller, the output power has an oscillating and unstable response. With the application of HHO-PID, the output power fluctuations are reduced to zero in less than 10 seconds.Figure 6
Output power before applying the controller.Figure 7
Output power after applying the controller.Figure 8
Nyquist diagram of the load angle.To evaluate the stability of the proposed controller, the Nyquist and Bode diagrams of the load angle are plotted in Figures8 and 9. As shown in Figure 8, this diagram is not circumvented, −1. According to the Nyquist stability criterion, since the open-loop transmission function does not have an instability pole, the load angle output is stable. Using Figure 9, since the gain margin and phase margin are positive, the HHO-PID controller will stabilize the load angle. Figures 8 and 9 clearly show the proposed controller performance in the load angle stability of the generator.Figure 9
Bode diagram of the load angle.Figure 10
Convergence curve of the proposed HHO-PID algorithm.


## 5.2. Proposed Controller Performance

As can be seen in Figure10, after 100 repetitions of the cost function, it has reached 27.1, which confirms the convergence of the proposed method. To evaluate the performance of the proposed controller, the step responses of the HHO-PID controller with fuzzy-PID (FPID) and sliding mode controllers (SMC) have been compared. Therefore, by applying a unit step as Ug input to the synchronous generator system, its load angle characteristic is obtained according to Figure 11(a). As shown in this figure, the characteristic of the load angle is wildly oscillating. Its overshoot and settling time are high, and, after 3000 seconds, it still oscillates, which is not desirable. As shown in Figure 11(e), the step response has a much better performance in its overshoot, settling time, fluctuations, and rise time.Figure 11
Load angle characteristic for a step input. (a) Without the controller, (b) with SMC in control level mode greater than one (SMCU), (c) with fuzzy-PID controller, (d) with SMC in control level mode less than one (SMCD), and (e) with HHO-PID controller.
(a)(b)(c)(d)(e)Also, the performances of the sliding mode controller (SMC) and the fuzzy-PID controller (FPID) are examined. By selecting the control level less than 1 in the Simulink model, the load angle output is obtained as in Figure11(d), which is not a suitable output. Figure 11(b) shows the load angle of a synchronous generator rotor with an SMC in a situation where the control level is greater than one. In this case, it can be seen that this controller does not work properly. Mamdani fuzzy controller with two inputs and three outputs as well as 49 rules is programmed. As shown in Figure 11(c), the load angle is completely unstable and therefore indicates that the FPID controller cannot be an excellent option to control the load angle of the synchronous generator rotor.In SMC and FPID controller, due to the fact that the power plant generator model uses mechanical parameters and the system order is 14, these controllers have not performed well and have diverged. In [28], the load angle of the 7th-order synchronous generator system is controlled by a PID controller using a genetic algorithm. A comparison of the performances of these GA-PID [28], Modified GA-PID, SMCU, SMCD, and FPID methods is presented in Table 4. The results show that the proposed algorithm and system are reduced by 89% at rise time, 72.6% at settling time, and 1.59% at overshoot and have better performance than GA-PID in all step response parameters. In the following, the optimization algorithm used is compared with the modified genetic algorithm. As can be seen in Table 4, the rise time and settling time of the proposed algorithm have been reduced by 36.6% and 74.2% compared to the genetic algorithm, respectively.Table 4
Comparison of the proposed controller performance with GA-PID [28], Modified GA-PID, SMCU, SMCD, and FPID.
Rise time (sec)Settling time (sec)Overshoot (%)Peak time (sec)HHO-PID0.436913.700854.12641.2101GA-PID [28]450558Modified GA-PID0.689753.187852.66781.9944SMCU1500——2000SMCD1022004820FPID2400——2440


## 6. Conclusion

In this paper, the HHO-PID controller is analyzed to improve the load angle stability and the characteristic of active power in the synchronous generator. The PID controller coefficients are optimized according to the objective function defined by the HHO algorithm and applied to the 14-order system of synchronous generators. Compared with other articles and methods, it is shown that the proposed algorithm has low rise time, low settling time, low peak time, and good overshoot in controlling the load angle and output power of the generator. It should be noted that, in this paper, the parameters of the HHO algorithm similar to other existing papers are calculated by trial and error. In future studies, researchers can work on optimizing the parameters of the HHO optimization algorithm.



---

*Source: 7746062-2022-03-14.xml*