# Multiobjective Optimal Control of FOPID Controller for Hydraulic Turbine Governing Systems Based on Reinforced Multiobjective Harris Hawks Optimization Coupling with Hybrid Strategies

Authors: Wenlong Fu; QiPeng Lu
Journal: Complexity
 (2020)

---

## Abstract

The controlling parameter tuning of the hydraulic turbine governing system (HTGS) is always deduced under single operating condition and is not suitable for the changeable operating conditions of the hydraulic turbine. For this purpose, multiobjective optimization problem of fractional order PID (FOPID) controller for HTGS is constructed through the consideration of no-load disturbance and on-load disturbance operation conditions, where the performance indicators of integral time absolute error (ITAE) under both operation conditions are employed as the objective functions. To achieve the optimum, the multiobjective version of newly proposed Harris hawks optimization (MOHHO) is established to solve the optimization issue. Additionally, hybrid strategies which include Latin hypercube sampling initialization, modified differential evolution operator, and mutation operator are coupled into MOHHO (HMOHHO) to promote the global searching capability. Simultaneously, the linear model of rabbit energy within MOHHO is replaced with a nonlinear one to further enhance the searching capacity. Subsequently, the effectiveness and superiority of the proposed HMOHHO are verified by several multiobjective UF and ZDT test problems. Finally, the practical application and contrastive analysis ascertain that the constructed multiobjective problem of FOPID controller is suitable for HTGS under changeable operating conditions, and the proposed HMOHHO is effective in solving the issue.


---

## Body

## 1. Introduction

The Energy Internet [1] is called the core technology of “the third industrial revolution,” which aims to replace fossil energy with renewable energy [2] and has attracted sufficient attention worldwide. The Energy Internet heightens the proportion of renewable energy gradually in primary energy production and consumption, thus to establish a sustainable energy supply system. Generally, renewable energy contains regular energy and new energy, where new energy has the characteristics of intermittence, randomness, and poor adjustment ability. By contrast, the main regular energy-hydropower energy possesses the features of low production cost, flexible operation, and no pollution [3]. Besides, with the increasing complexity of energy structure, hydropower undertakes more and more tasks of peak load regulation [4] and frequency modulation [5]. In this situation, the operation safety, reliability, and health management [6] of hydropower generator [7–9] which is the key device during the hydropower energy conversion process are of particular importance.In the process of frequent working condition transformation, hydraulic turbine governing system (HTGS) acts as the main role for ensuring the effective operation of hydropower generator. Nevertheless, there always exist strong nonlinear characteristics within HTGS, which would affect the dynamic performance greatly [4, 10]. Hence, favourable control strategy is very critical to promote the performance of HTGS. In previous studies, there are various control methods which have been used to solve control problems in HTGS, such as sliding mode control [11–14], fault tolerance control [15], and predictive control [16]. In the controlling field of HTGS, the most common control strategy is traditional PID control strategy which is uncomplicated and is simple to implement [17]. For example, Jiang et al. [18] proposed an optimized PID controller for hydroturbine governing system. Khodabakhshian and Hooshmand [19] proposed a new robust PID controller for automatic generation control of hydroturbine power systems. However, the adjusting ability of traditional PID controller is limited, which would result in poor performance of HTGS under multiple operating conditions. By contrast, the FOPID control strategy can achieve adjustment flexibility for hydropower unit [20], which has been studied by researchers. For example, Wu et al. [21] proposed a fuzzy fractional order PID (FFOPID) controller through the combination of FOPID and fuzzy logic controller. The experimental results showed that the optimized design of FFOPID controller has better control quality than PID, FOPID, and fuzzy PID (FPID) controller. Xu et al. [22] proposed a robust nonfragile fractional order PID controller to solve pumped turbine unit frequency oscillation along the “S” shape curve under unload operation condition. The simulation experiment verified that the proposed FOPID controller based on the bacterial foraging algorithm had better robustness and practicability than PID controller under no-load operation condition. Li et al. [23] designed a FOPID controller for a pumped storage unit, and results showed that the FOPID controller optimized by gravitational search algorithm combined with Cauchy and Gaussian mutation obtained obvious advantage over other PID controllers.Apart from the superiority of FOPID controller over PID controller, the parameter tuning of controller is an extremely important step for employment of FOPID controller in HTGS. With the sustainable development of intelligence algorithms in recent years, intelligent algorithms have been widely applied to optimize the control parameters of HTGS by searching the defined available space. The commonly used algorithms include the following: genetic algorithm (GA) [24], particle swarm optimization (PSO) [25, 26], gravitational search algorithm (GSA) [27, 28], ant lion optimization (ALO) [29], sine cosine algorithm (SCA) [30, 31], and grey wolf optimization (GWO) [32–34]. Meanwhile, error indicator is always employed to measure the performance when optimizing control parameters with intelligent algorithms. Generally, the parameter selection considering multiple error indexes can make up for the shortcomings of single error index. In particular, with regard to two commonly used indexes, integral squared error (ISE) and integral time squared error (ITSE) which are opposed, system applying ISE index has a response with small overshoot percentage but long settling time, while system applying ITSE index has a response with shorter settling time but without stability margin. In other words, the system designer can set the error index weight percentage based on specific system requirements. For example, Chen et al. [35] proposed chaotic nondominated sorting genetic algorithm II (NSGAII) to optimize a multiobjective optimization problem of FOPID controller’s parameter tuning, whose objective functions are composed of ISE and ITSE. Piraisoodi et al. [36] proposed a multiobjective robust fuzzy FOPID controller designed for nonlinear HTGS by using NSGAII. The results showed that the proposed controller had better fitness value and time domain specifications than PID and FOPID controllers, as well as satisfying the conflicting objectives including less settling time and minimum damped oscillations. However, their control parameters are always optimized by intelligent algorithms under a single operating condition, in which case the optimization results may not be adaptive for the changeable operating conditions of hydropower unit.To promote the feasibility of HTGS under changeable working conditions, it is rather necessary and important to consider the optimal control parameters of HTGS under multiple operating conditions [37]. For example, Zhang et al. [38] proposed an improved NSGAIII algorithm to solve multiobjective FOPID controller optimization problem for pumped turbine governing system (PTGS) under multiworking conditions. Xia et al. [39] proposed a multiobjective PID controller for HTGS based on an improved MOGWO algorithm under multiworking conditions. In this paper, the integral time absolute error (ITAE) indexes of HTGS operating at changeable working conditions are constructed to deduce the optimal control parameters of HTGS under different working conditions. Essentially, the control parameters optimization of HTGS under different working conditions could be summarized as a multiobjective optimization issue, which is expected to be solved by multiobjective optimization algorithms. Some representative multiobjective optimization algorithms include the following: nondominated sorting genetic algorithm III (NSGA-III) [40], multiobjective particle swarm optimization (MOPSO) [41], and multiobjective grey wolf optimizer (MOGWO) [42]. Among the above algorithms, MOPSO and MOGWO are the multiobjective versions of the corresponding single objective algorithms. Inspired by this condition, a multiobjective version of Harris hawks optimization (MOHHO) is structured based on single objective HHO, which is recently proposed by Heidari et al. [43] in 2019 and whose superiority has been ascertained. Furthermore, the multiobjective HHO is reinforced with hybrid strategies (HMOHHO) including Latin hypercube sampling initialization, modified differential evolution operator, and mutation operator, to deduce the optimal control parameters of FOPID controller for HTGS under multiworking conditions. Simultaneously, the linear model of rabbit energy within MOHHO is replaced with a nonlinear one to further enhance the searching capacity.The remainder of this paper is organized as follows: Section2 briefly introduces the basic concepts of fractional calculus and FOPID controller. In Section 3, HTGS model and its control issues are discussed. The proposed HMOHHO algorithm is constructed in Section 4. In Section 5, the proposed HMOHHO is compared with NSGAIII, MOPSO, MOGWO, and MOHHO and performance analysis is accomplished. In Section 6, comparison of control effects between PID and FOPID controllers is described. Then the Pareto-optimal sets with different algorithms are compared and analyzed, and the transient process responses with different control schemes are presented. Lastly, Section 7 summarizes the conclusions.


## 2. Fractional Calculus and the Fractional Order PID

### 2.1. Theory of Fractional Calculus

The research of fractional calculus [44] has been carried out since 1960s and widely expanded into the fields of science and technology. During the practical applications, there are always some complex issues which make it difficult to explain the model with integer order calculus, while fractional calculus shows more flexibility and availability. For this purpose, the theory of fractional calculus is introduced and then employed to promote the controller modeling for HTGS later. Here, a unified fractional calculus operator (t0Dtα) is introduced and defined as(1)t0Dtαft=dαdtαft,α>0,ft,α=0,∫t0tfτdτ−α,α<0,where α is limited to real number, t is an independent variable, and t0 is the lower boundary of t.The definition of fractional calculus was firstly proposed in 1868, which meant the real establishment in the field of fractional calculus. Compared with the definitions of Riemann–Liouville and Grünwald–Letnikov [45], the definition proposed by Caputo [46] is more suitable for conditions with nonzero initialization. The function of fractional derivative defined by Caputo is(2)Dtαt0cyt=1Γm−α∫t0tymτt−τ1+α−mdτ,where α∈ R, m ∈ Z, m − 1 < α ≤ m, γ is limited to real number, and Γ(·) is Euler Gamma function.As it can be seen from formula (2), Caputo’s definition requires that the m-th derivatives of the function are integrable [47]. The Laplace transform of fractional derivative defined by Caputo is expressed in(3)∫0∞e−stDγftdt=sγFs−∑k=0n−1sγ−k−1fk0,where n is the smallest integer, γ donates the order of fractional derivative, and n − 1 < γ ≤ n.On account of unappealing outcomes for frequency response fitting with design method of filter based on continued fraction, the filter of fractional order operator proposed by Oustaloup can choose the frequency bands and order, which can approximate the fractional calculus operator with integer order transfer function. However, approximation effect of Oustaloup filter at selected frequency bands boundary is unsatisfactory. This article uses an improved filter form Oustaloup filter [47], whose mathematical model is(4)sγ≈dwhbγds2+bωhsd1−γs2+vwhs+dγ∏k=1Ns+ωk′s+ωk.In formula (4), the zero point, pole point, and gain can be calculated, respectively, as follows:(5)ωu=ωhωb,ωk′=ωbωu2k−1−γ/N,ωk=ωbωu2k−1−γ/N,K=ωhγ,where (ωb, ωh) denotes frequency band, γ ∈ (0, 1). In general, the weighting parameter is selected as b = 10 and d = 9. The filter order N is set at 13 and the frequency band is set (10−5, 103) in this paper.


### 2.2. Basic Concepts of FOPID Controller

The parametersλ and μ of FOPID controller can be set any real number between 0 and 2, which is a generalized form of traditional integer order PID controller. Compared with the traditional integer order PID controller, the parameter tuning of FOPID controller increases the complexity of algorithmic calculation process to a certain extent due to the additional two parameters. However, FOPID controller plays a great important role in improving the flexibility, robustness, and overall control effect of system.The mathematical models of FOPID controller including the time-domain model and frequency-domain transfer function model are(6)ut=KPet+KiD−λet+KdDμet,Cs=Kp+Kis−λ+Kdsμ,λ,μ>0,where λ = 1 and μ = 1; it is the traditional PID controller model.


## 2.1. Theory of Fractional Calculus

The research of fractional calculus [44] has been carried out since 1960s and widely expanded into the fields of science and technology. During the practical applications, there are always some complex issues which make it difficult to explain the model with integer order calculus, while fractional calculus shows more flexibility and availability. For this purpose, the theory of fractional calculus is introduced and then employed to promote the controller modeling for HTGS later. Here, a unified fractional calculus operator (t0Dtα) is introduced and defined as(1)t0Dtαft=dαdtαft,α>0,ft,α=0,∫t0tfτdτ−α,α<0,where α is limited to real number, t is an independent variable, and t0 is the lower boundary of t.The definition of fractional calculus was firstly proposed in 1868, which meant the real establishment in the field of fractional calculus. Compared with the definitions of Riemann–Liouville and Grünwald–Letnikov [45], the definition proposed by Caputo [46] is more suitable for conditions with nonzero initialization. The function of fractional derivative defined by Caputo is(2)Dtαt0cyt=1Γm−α∫t0tymτt−τ1+α−mdτ,where α∈ R, m ∈ Z, m − 1 < α ≤ m, γ is limited to real number, and Γ(·) is Euler Gamma function.As it can be seen from formula (2), Caputo’s definition requires that the m-th derivatives of the function are integrable [47]. The Laplace transform of fractional derivative defined by Caputo is expressed in(3)∫0∞e−stDγftdt=sγFs−∑k=0n−1sγ−k−1fk0,where n is the smallest integer, γ donates the order of fractional derivative, and n − 1 < γ ≤ n.On account of unappealing outcomes for frequency response fitting with design method of filter based on continued fraction, the filter of fractional order operator proposed by Oustaloup can choose the frequency bands and order, which can approximate the fractional calculus operator with integer order transfer function. However, approximation effect of Oustaloup filter at selected frequency bands boundary is unsatisfactory. This article uses an improved filter form Oustaloup filter [47], whose mathematical model is(4)sγ≈dwhbγds2+bωhsd1−γs2+vwhs+dγ∏k=1Ns+ωk′s+ωk.In formula (4), the zero point, pole point, and gain can be calculated, respectively, as follows:(5)ωu=ωhωb,ωk′=ωbωu2k−1−γ/N,ωk=ωbωu2k−1−γ/N,K=ωhγ,where (ωb, ωh) denotes frequency band, γ ∈ (0, 1). In general, the weighting parameter is selected as b = 10 and d = 9. The filter order N is set at 13 and the frequency band is set (10−5, 103) in this paper.


## 2.2. Basic Concepts of FOPID Controller

The parametersλ and μ of FOPID controller can be set any real number between 0 and 2, which is a generalized form of traditional integer order PID controller. Compared with the traditional integer order PID controller, the parameter tuning of FOPID controller increases the complexity of algorithmic calculation process to a certain extent due to the additional two parameters. However, FOPID controller plays a great important role in improving the flexibility, robustness, and overall control effect of system.The mathematical models of FOPID controller including the time-domain model and frequency-domain transfer function model are(6)ut=KPet+KiD−λet+KdDμet,Cs=Kp+Kis−λ+Kdsμ,λ,μ>0,where λ = 1 and μ = 1; it is the traditional PID controller model.


## 3. Multiobjective Optimization Framework for FOPID and PID Controllers in HTGS

### 3.1. HTGS Model

HTGS system, including governor, hydraulic servo system, hydraulic turbine, penstock, and generator [48], is influenced by hydraulic, mechanical, and electrical factors, making its response behavior complicated. The structure diagram of HTGS is shown in Figure 1. As far as system modeling is concerned, the simulation system expressed roundly by various factors plays a significant role in promoting the reliability of simulation results. Nevertheless, in the modeling of actual complex systems, there always exist many influencing factors which are difficult to be considered totally, which means that the least important factors are neglected in most cases.Figure 1
The structure diagram of HTGS.During the transition process of hydropower station, the change of flowing water in the penstock would induce water hammer effect. When the length of penstock is less than 800 meters, it is considered that the elasticity of water body and penstock wall has little effect on the water hammer, amounting to that water hammer pressure spreading to the entire penstock is accomplished instantaneously. The rigid water hammer model of water diversion system is(7)Ghs=hsqs=−Tws,where Tw is the water flow inertia time constant, which is an important parameter in the water diversion system.Torque and flow rate of the Francis hydraulic turbine are related to the guide vane opening, rotational speed, and water head. The Francis hydraulic turbine model in steady state can be expressed as(8)mt=mty,h,ω,q=qy,h,ω,where mt, q, ω, y, and h, respectively, represent relative torque deviation relative value, flow deviation relative value, rotational speed deviation relative value, guide vane opening deviation relative value, and head deviation relative value.The Taylor expansion of equation (8) is simplified by omitting the second order and higher order differential components [49]. Thus, the following formula can be obtained:(9)mt=exω+eyy+ehh,q=eqxω+eqyy+eqhh,where ex denotes the first order partial derivative value of torque in relation to speed of hydraulic turbine, ey denotes first order partial derivative value of torque in relation to wicket gate, eh denotes the first order partial derivative value of torque with respect to water head, eqx denotes the first order partial derivative value of flow rate in relation to speed of hydraulic turbine, eqy denotes the first order partial derivative value of flow rate in relation to wicket gate, and eqh denotes the first order partial derivative value of flow rate in relation to water head.The generator is also a complicated subsystem, which is divided into different models according to the order of differential equation. In this paper, the first order generator model is introduced and researched:(10)Ggs=xsmts−mgs=1Tas+eg,where x donates the frequency of generator, mg represents load torque relative deviation, Ta represents generator mechanical time, and eg represents generator load self-regulation parameters.Dead zones and amplitude limit units have been considered in hydraulic servo subsystem. The transfer function of the hydraulic actuator can be expressed as(11)Gys=1Tys+1.where Ty is the major relay connecter response time.According to the mathematical model described above, the system block diagram of hydraulic turbine unit can be modeled as shown in Figure2.Figure 2
The system block diagram of hydraulic turbine unit.The governor is the key equipment for automation of hydropower station. In other words, the control strategy directly affects the safety and stability of hydropower station and unit. The block diagrams of PID and FOPID controller (the amplitude limit unit is ignored) are shown in Figure3.Figure 1
Controllers: (a) PID and (b) FOPID.
(a)(b)


### 3.2. Multiobjective Optimal Control of HTGS and Problem Description

During the process of generating electricity in hydropower station, speed control of HTGS is an extremely important link of the unit control and automation, the quality of which is affected by the parameter tuning of controller. In essence, the parameter selection of controller is an optimization problem to solve extreme values. According to the system requirements, objective function is set with the error performance indicator firstly. Then, a certain method is used to optimize the parameters of controller by objective function. The most commonly used performance indicators [35, 50] are ISE, integral absolute error (IAE), ITSE, ITAE, integral squared time squared error (ISTSE), and integral squared time absolute error (ISTAE).The dynamic responding performance of HTGS is vital for the robustness and stability of power system, which is influenced by load fluctuation and severe frequency interference. Thus, it is expected that the performance objective functions can achieve accurate and robust tracking control. In this paper, ITAE which is one of the most widely used indicators with stable regulation and small overshoot is chosen as the objective function under multiple operating conditions to obtain better transient dynamic performance of system. The ITAE index [35] is defined as follows:(12)ITAE=∫0Ttetdt,where e(t) denotes the relative deviation of the rotational speed within HTGS.Considering two classic operating conditions, single objective HHO is applied to optimize parameters of FOPID controller under 4% step disturbance condition (no-load operation condition) with the transfer parameters in Table1 and then HTGS is controlled with the obtained parameters under 4% load shedding condition (on-load operation condition). The parameters are tuned as Kp = 3.5831, Ki = 1.1739, Kd = 4.5975, λ = 1.0009, and μ = 0.0010. Besides, the control effects of tuned parameters under two operating conditions are shown in Figure 4, from which it can be seen that the parameters of FOPID controller obtained under no-load disturbance are not suitable for the on-load disturbance. In other words, the parameters of controller obtained under single operating condition are not desirable for the control of hydraulic turbine unit under variable operating conditions.Table 1
Transferring parameters in hydroturbine and generator under two working conditions.
Running conditionexeyeheqxeqyeqhegTaTwTyK1vOn-load−0.3121.1120.4680.356−0.4020.3010.8980.560.10.28No-load−1.0041.3411.2930.3691.0740.2970.8980.560.10.28Figure 4
Control effects of the same parameter under two operating conditions. (a) No-load operating condition; (b) on-load operating condition.
(a)(b)The main reason of conclusion above is that no-load disturbance and on-load disturbance operation condition are extreme operating conditions in operating process of hydraulic turbine unit. Therefore, optimal control of HTGS is essentially a multiobjective optimization problem. Different operating conditions are considered in this paper for optimal control. Referring to [38], the objective functions of HTGS based on ITAE under two conditions are described as follows:(13)minf1=ITAE1=f1Kp,Ki,Kd,λ,μf2=ITAE2=f2Kp,Ki,Kd,λ,μsubject toKpmin≤Kp≤KpmaxKimin≤Ki≤KimaxKdmin≤Kd≤Kdmaxλmin≤λ≤λmaxμmin≤μ≤μmax,where f1(·) and f2(·) are functions of Kp, Ki, Kd, λ, and μ under on-load and no-load disturbance running conditions; the lower and upper bounds of Kp, Ki, Kd, λ, and μ are Kpmin = 0, Kpmax = 10, Kimin = 0, Kimax = 10, Kdmin = 0, Kdmax = 10, λmin = 0, λmax = 2, μmin = 0, and μmax = 2.


## 3.1. HTGS Model

HTGS system, including governor, hydraulic servo system, hydraulic turbine, penstock, and generator [48], is influenced by hydraulic, mechanical, and electrical factors, making its response behavior complicated. The structure diagram of HTGS is shown in Figure 1. As far as system modeling is concerned, the simulation system expressed roundly by various factors plays a significant role in promoting the reliability of simulation results. Nevertheless, in the modeling of actual complex systems, there always exist many influencing factors which are difficult to be considered totally, which means that the least important factors are neglected in most cases.Figure 1
The structure diagram of HTGS.During the transition process of hydropower station, the change of flowing water in the penstock would induce water hammer effect. When the length of penstock is less than 800 meters, it is considered that the elasticity of water body and penstock wall has little effect on the water hammer, amounting to that water hammer pressure spreading to the entire penstock is accomplished instantaneously. The rigid water hammer model of water diversion system is(7)Ghs=hsqs=−Tws,where Tw is the water flow inertia time constant, which is an important parameter in the water diversion system.Torque and flow rate of the Francis hydraulic turbine are related to the guide vane opening, rotational speed, and water head. The Francis hydraulic turbine model in steady state can be expressed as(8)mt=mty,h,ω,q=qy,h,ω,where mt, q, ω, y, and h, respectively, represent relative torque deviation relative value, flow deviation relative value, rotational speed deviation relative value, guide vane opening deviation relative value, and head deviation relative value.The Taylor expansion of equation (8) is simplified by omitting the second order and higher order differential components [49]. Thus, the following formula can be obtained:(9)mt=exω+eyy+ehh,q=eqxω+eqyy+eqhh,where ex denotes the first order partial derivative value of torque in relation to speed of hydraulic turbine, ey denotes first order partial derivative value of torque in relation to wicket gate, eh denotes the first order partial derivative value of torque with respect to water head, eqx denotes the first order partial derivative value of flow rate in relation to speed of hydraulic turbine, eqy denotes the first order partial derivative value of flow rate in relation to wicket gate, and eqh denotes the first order partial derivative value of flow rate in relation to water head.The generator is also a complicated subsystem, which is divided into different models according to the order of differential equation. In this paper, the first order generator model is introduced and researched:(10)Ggs=xsmts−mgs=1Tas+eg,where x donates the frequency of generator, mg represents load torque relative deviation, Ta represents generator mechanical time, and eg represents generator load self-regulation parameters.Dead zones and amplitude limit units have been considered in hydraulic servo subsystem. The transfer function of the hydraulic actuator can be expressed as(11)Gys=1Tys+1.where Ty is the major relay connecter response time.According to the mathematical model described above, the system block diagram of hydraulic turbine unit can be modeled as shown in Figure2.Figure 2
The system block diagram of hydraulic turbine unit.The governor is the key equipment for automation of hydropower station. In other words, the control strategy directly affects the safety and stability of hydropower station and unit. The block diagrams of PID and FOPID controller (the amplitude limit unit is ignored) are shown in Figure3.Figure 1
Controllers: (a) PID and (b) FOPID.
(a)(b)


## 3.2. Multiobjective Optimal Control of HTGS and Problem Description

During the process of generating electricity in hydropower station, speed control of HTGS is an extremely important link of the unit control and automation, the quality of which is affected by the parameter tuning of controller. In essence, the parameter selection of controller is an optimization problem to solve extreme values. According to the system requirements, objective function is set with the error performance indicator firstly. Then, a certain method is used to optimize the parameters of controller by objective function. The most commonly used performance indicators [35, 50] are ISE, integral absolute error (IAE), ITSE, ITAE, integral squared time squared error (ISTSE), and integral squared time absolute error (ISTAE).The dynamic responding performance of HTGS is vital for the robustness and stability of power system, which is influenced by load fluctuation and severe frequency interference. Thus, it is expected that the performance objective functions can achieve accurate and robust tracking control. In this paper, ITAE which is one of the most widely used indicators with stable regulation and small overshoot is chosen as the objective function under multiple operating conditions to obtain better transient dynamic performance of system. The ITAE index [35] is defined as follows:(12)ITAE=∫0Ttetdt,where e(t) denotes the relative deviation of the rotational speed within HTGS.Considering two classic operating conditions, single objective HHO is applied to optimize parameters of FOPID controller under 4% step disturbance condition (no-load operation condition) with the transfer parameters in Table1 and then HTGS is controlled with the obtained parameters under 4% load shedding condition (on-load operation condition). The parameters are tuned as Kp = 3.5831, Ki = 1.1739, Kd = 4.5975, λ = 1.0009, and μ = 0.0010. Besides, the control effects of tuned parameters under two operating conditions are shown in Figure 4, from which it can be seen that the parameters of FOPID controller obtained under no-load disturbance are not suitable for the on-load disturbance. In other words, the parameters of controller obtained under single operating condition are not desirable for the control of hydraulic turbine unit under variable operating conditions.Table 1
Transferring parameters in hydroturbine and generator under two working conditions.
Running conditionexeyeheqxeqyeqhegTaTwTyK1vOn-load−0.3121.1120.4680.356−0.4020.3010.8980.560.10.28No-load−1.0041.3411.2930.3691.0740.2970.8980.560.10.28Figure 4
Control effects of the same parameter under two operating conditions. (a) No-load operating condition; (b) on-load operating condition.
(a)(b)The main reason of conclusion above is that no-load disturbance and on-load disturbance operation condition are extreme operating conditions in operating process of hydraulic turbine unit. Therefore, optimal control of HTGS is essentially a multiobjective optimization problem. Different operating conditions are considered in this paper for optimal control. Referring to [38], the objective functions of HTGS based on ITAE under two conditions are described as follows:(13)minf1=ITAE1=f1Kp,Ki,Kd,λ,μf2=ITAE2=f2Kp,Ki,Kd,λ,μsubject toKpmin≤Kp≤KpmaxKimin≤Ki≤KimaxKdmin≤Kd≤Kdmaxλmin≤λ≤λmaxμmin≤μ≤μmax,where f1(·) and f2(·) are functions of Kp, Ki, Kd, λ, and μ under on-load and no-load disturbance running conditions; the lower and upper bounds of Kp, Ki, Kd, λ, and μ are Kpmin = 0, Kpmax = 10, Kimin = 0, Kimax = 10, Kdmin = 0, Kdmax = 10, λmin = 0, λmax = 2, μmin = 0, and μmax = 2.


## 4. Reinforced Multiobjective Harris Hawks Optimization with Hybrid Strategies

### 4.1. Harris Hawks Optimization

Mainly inspired by the chasing style and cooperative behavior of Harris hawks, Heidari and his team developed HHO algorithm which has two main stages: exploration stage and exploitation stage [43]. During exploration stage, Harris’s hawks inhabit in random places, waiting to detect prey based on two strategies. The mathematical model is as follows:(14)Xt+1=Xrandt−r1Xrandt−2r2Xt,q≥0.5,Xrabbitt−Xmt−r3LB+r4UB−LB,q<0.5,where X(t + 1) represents the position of hawks in t + 1th iteration, Xrand(t) donates the position of hawks chosen randomly in the current population, r1, r2, r3, r4, and q are all random numbers in scope of (0, 1), Xrabbit(t) is the position vector of rabbit, Xt represents the current position of hawk, which is updated during each iteration, LB and UB denote the upper and lower limits of variable, and Xm(t) represents the average position of hawks in the current population, which can be calculated as(15)Xmt=1N∑i=1NXit,where Xi(t) denotes the location of each hawk in t-th iteration and N represents the amount of hawks.During the transition phase from exploration to exploitation, the energy of rabbit is modeled as(16)E=2E01−tT,where E is the evasion energy of rabbit, T denotes the maximum iterations, and E0 is a value in scope of (−1, 1), indicating the initial energy of each step.WhenE > 1, optimization process of Harris hawks is focused mainly in the exploration stage; otherwise, it turns to the exploitation stage which includes the stages of soft besiege, hard besiege, soft besiege with progressive rapid dives, and hard besiege with progressive rapid dives.In soft besiege stager≥0.5andE≥0.5, the behavioral model can be constructed as follows:(17)Xt+1=ΔXt−EJXrabbitt−Xt,ΔXt=Xrabbitt−Xt,where ΔX(t) is the difference value between the position of rabbit and the current position in iteration t, J = 2(1 − r5) denotes the random energy of rabbit during the escaping procedure, and r5 is the random number in scope of (0, 1). The changed value of J is realized to simulate the trait of rabbit’s escaping movement.In hard besiege stager≥0.5andE≤0.5, the current position is updated by the following formula:(18)Xt+1=Xrabbitt−E△Xt.In soft besiege with progressive rapid dives stager<0.5andE≥0.5, in order to carry out advanced soft besiege, hawks can assess their next action according to the following rule:(19)Y=Xrabbitt−EJXrabbitt−Xt.The next step will also use the following rule for diving based on levy flight (LF) mode:(20)Z=Y+S×LFD,where D represents the dimension of decision vector, S represents a 1 × D random vector in scope of (−1, 1), and LF shows the levy flight function formulated as(21)LFx=0.01×u×σv1/β,σ=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β,β=1.5,where u and v are both random values in scope of (0, 1).Therefore, in order to get a better position vector, the last strategy of soft besiege with progressive rapid dives stage can be expressed as(22)Xt+1=Y,ifFY<FZ,Z,ifFZ<FY.In hard besiege with progressive rapid dives stager<0.5andE<0.5, the following rules are executed in hard besiege condition:(23)Xt+1=Y,ifFY<FZ,Z,ifFZ<FY,where Y and Z can be deduced by the following novel rules:(24)Y=Xrabbitt−EJXrabbitt−Xmt,Z=Y+S×LFD.


### 4.2. Reinforced MOHHO with Hybrid Strategies

#### 4.2.1. MOHHO

Multiobjective Harris hawks optimization algorithm (MOHHO) based on HHO is expected to solve multiobjective optimization problem. The steps of MOHHO algorithm are as follows:Step 1: preset the population sizeN, maximum iterations T, and archive size S at the beginning as well as initializing the population randomly within the given ranges of variablesStep 2: calculate the objective valuesf1(·) and f2(·) for each hawk X(t)Step 3: select the nondominated Pareto-optimal solutions in the population at the current iteration; thus, the position of rabbitXrabbit(t) is acquired by leader selection mechanism, whose role is taking advantage of crowding distance to select a solution through roulette wheel method from a less populated area of the archiveStep 4: update each hawkX(t) according to equations (14)–(24)Step 5: calculate the new objective value for each hawk and find the nondominated solutions; thus better solutions will be recorded in the archiveStep 7: when the archive is full, the crowded area of the archive is deleted by the roulette wheel method for adding new solutions to the archiveStep 8: output the archive solutions


#### 4.2.2. Reinforced MOHHO

Previous literatures [51–54] have shown that heuristic search algorithms generally suffer from the problem of trapping in local optimum easily, leading to the loss of solution diversity. Similarly, the original iterative process of MOHHO is not sufficient to maintain lateral diversity and achieve a Pareto front of good astringency and high diversity. In order to promote the searching capacity of MOHHO, the linear model of rabbit energy is replaced with a nonlinear one, as shown in(25)E=2E01−tT6.In addition to the improvement above, hybrid strategies are merged into MOHHO to jump out local optimum and search for more no-domain solutions. Firstly, population initialization is conducted through Latin hypercube sampling, which is the latest development in sampling technology [55]. The steps of Latin hypercube sampling initialization are as follows: (1) divide each dimension into N intervals that do not overlap each other; (2) randomly select a point within each interval for each dimension; (3) combine them into a vector. The mathematical model can be described by the following formula:(26)xd_min=xd0<xd1<xd2<⋯<xdj<⋯<xdN=xd_max,where xd donates the d-th dimension of each hawk vector and N donates population size.The modified differential evolution operator of the proposed hybrid strategies is shown in(27)Xt+1=Xrand1+F∗Xrand2−Xrand3,where Xrand(1), Xrand(2), and Xrand(3) are the positions of hawks chosen randomly in the current population, F is a 1 × dim vector whose value intervals for all dimensions are 0.3 and 0.8, and dim is the dimension number of each hawk.The mutation operator mathematical model of the proposed hybrid strategies can be expressed as(28)d=m∗UBj−LBj,j∈rand1,2,…,dim,Xt+1=Xj=Xj+d,where m is a random value in scope of (−0.5, 0.5).The specific pseudocode of the proposed HMOHHO algorithm is depicted in Algorithm1.Algorithm 1: Pseudocode of HMOHHO algorithm.
Inputs: Population size N, maximum iterations T and archive size SInitialize: Latin hypercube sampling initialize population Xii=1,2,…,NCalculate f1(·) and f2(·)Find the non-dominated solutionswhilet < Tifthe archive is fullCrowded area of the archive is deleted by roulette wheel methodAdd new solutions to the archiveend ifSelect Xrabbit(t) in the archiveforeach hawkifabs (E) < 1X(t) is updated by (14)elseX(t) is updated by (15)–(24)end ifXr1 is updated by (27)iff(Xr1) dominate f(X(t))X(t) = Xr1end ifXr2 is updated by (28)iff(Xr2) dominate f(X(t))X(t) = Xr2end ifend forCalculate f1(·) and f2(·)Find the non-dominated solutionst = t + 1end whilereturnarchive


### 4.3. Computational Complexity Analysis

The complexity analysis of each iteration in the proposed HMOHHO is as follows: assume that individual number isN and the dimension of each hawk is dim, and the number of individuals in the current archive set is assumed to be A. In HMOHHO, the rabbit position in archive set is selected by leader selection mechanism, during which process the computational complexity can be denoted by O(A). Then, computational complexity of location updating can be denoted by O(N∗ dim). Additionally, the computational complexity of the proposed hybrid strategies can be denoted by O(2 ∗N∗ dim) and the computational complexity of updating archive set can be denoted by O(A ∗N∗ dim). Next, the computational complexity of grouping and sorting archive sets is O(A2). On the whole, the computational complexity of each iteration in HMOHHO is O(A2) +O(A∗N∗ dim).


## 4.1. Harris Hawks Optimization

Mainly inspired by the chasing style and cooperative behavior of Harris hawks, Heidari and his team developed HHO algorithm which has two main stages: exploration stage and exploitation stage [43]. During exploration stage, Harris’s hawks inhabit in random places, waiting to detect prey based on two strategies. The mathematical model is as follows:(14)Xt+1=Xrandt−r1Xrandt−2r2Xt,q≥0.5,Xrabbitt−Xmt−r3LB+r4UB−LB,q<0.5,where X(t + 1) represents the position of hawks in t + 1th iteration, Xrand(t) donates the position of hawks chosen randomly in the current population, r1, r2, r3, r4, and q are all random numbers in scope of (0, 1), Xrabbit(t) is the position vector of rabbit, Xt represents the current position of hawk, which is updated during each iteration, LB and UB denote the upper and lower limits of variable, and Xm(t) represents the average position of hawks in the current population, which can be calculated as(15)Xmt=1N∑i=1NXit,where Xi(t) denotes the location of each hawk in t-th iteration and N represents the amount of hawks.During the transition phase from exploration to exploitation, the energy of rabbit is modeled as(16)E=2E01−tT,where E is the evasion energy of rabbit, T denotes the maximum iterations, and E0 is a value in scope of (−1, 1), indicating the initial energy of each step.WhenE > 1, optimization process of Harris hawks is focused mainly in the exploration stage; otherwise, it turns to the exploitation stage which includes the stages of soft besiege, hard besiege, soft besiege with progressive rapid dives, and hard besiege with progressive rapid dives.In soft besiege stager≥0.5andE≥0.5, the behavioral model can be constructed as follows:(17)Xt+1=ΔXt−EJXrabbitt−Xt,ΔXt=Xrabbitt−Xt,where ΔX(t) is the difference value between the position of rabbit and the current position in iteration t, J = 2(1 − r5) denotes the random energy of rabbit during the escaping procedure, and r5 is the random number in scope of (0, 1). The changed value of J is realized to simulate the trait of rabbit’s escaping movement.In hard besiege stager≥0.5andE≤0.5, the current position is updated by the following formula:(18)Xt+1=Xrabbitt−E△Xt.In soft besiege with progressive rapid dives stager<0.5andE≥0.5, in order to carry out advanced soft besiege, hawks can assess their next action according to the following rule:(19)Y=Xrabbitt−EJXrabbitt−Xt.The next step will also use the following rule for diving based on levy flight (LF) mode:(20)Z=Y+S×LFD,where D represents the dimension of decision vector, S represents a 1 × D random vector in scope of (−1, 1), and LF shows the levy flight function formulated as(21)LFx=0.01×u×σv1/β,σ=Γ1+β×sinπβ/2Γ1+β/2×β×2β−1/21/β,β=1.5,where u and v are both random values in scope of (0, 1).Therefore, in order to get a better position vector, the last strategy of soft besiege with progressive rapid dives stage can be expressed as(22)Xt+1=Y,ifFY<FZ,Z,ifFZ<FY.In hard besiege with progressive rapid dives stager<0.5andE<0.5, the following rules are executed in hard besiege condition:(23)Xt+1=Y,ifFY<FZ,Z,ifFZ<FY,where Y and Z can be deduced by the following novel rules:(24)Y=Xrabbitt−EJXrabbitt−Xmt,Z=Y+S×LFD.


## 4.2. Reinforced MOHHO with Hybrid Strategies

### 4.2.1. MOHHO

Multiobjective Harris hawks optimization algorithm (MOHHO) based on HHO is expected to solve multiobjective optimization problem. The steps of MOHHO algorithm are as follows:Step 1: preset the population sizeN, maximum iterations T, and archive size S at the beginning as well as initializing the population randomly within the given ranges of variablesStep 2: calculate the objective valuesf1(·) and f2(·) for each hawk X(t)Step 3: select the nondominated Pareto-optimal solutions in the population at the current iteration; thus, the position of rabbitXrabbit(t) is acquired by leader selection mechanism, whose role is taking advantage of crowding distance to select a solution through roulette wheel method from a less populated area of the archiveStep 4: update each hawkX(t) according to equations (14)–(24)Step 5: calculate the new objective value for each hawk and find the nondominated solutions; thus better solutions will be recorded in the archiveStep 7: when the archive is full, the crowded area of the archive is deleted by the roulette wheel method for adding new solutions to the archiveStep 8: output the archive solutions


### 4.2.2. Reinforced MOHHO

Previous literatures [51–54] have shown that heuristic search algorithms generally suffer from the problem of trapping in local optimum easily, leading to the loss of solution diversity. Similarly, the original iterative process of MOHHO is not sufficient to maintain lateral diversity and achieve a Pareto front of good astringency and high diversity. In order to promote the searching capacity of MOHHO, the linear model of rabbit energy is replaced with a nonlinear one, as shown in(25)E=2E01−tT6.In addition to the improvement above, hybrid strategies are merged into MOHHO to jump out local optimum and search for more no-domain solutions. Firstly, population initialization is conducted through Latin hypercube sampling, which is the latest development in sampling technology [55]. The steps of Latin hypercube sampling initialization are as follows: (1) divide each dimension into N intervals that do not overlap each other; (2) randomly select a point within each interval for each dimension; (3) combine them into a vector. The mathematical model can be described by the following formula:(26)xd_min=xd0<xd1<xd2<⋯<xdj<⋯<xdN=xd_max,where xd donates the d-th dimension of each hawk vector and N donates population size.The modified differential evolution operator of the proposed hybrid strategies is shown in(27)Xt+1=Xrand1+F∗Xrand2−Xrand3,where Xrand(1), Xrand(2), and Xrand(3) are the positions of hawks chosen randomly in the current population, F is a 1 × dim vector whose value intervals for all dimensions are 0.3 and 0.8, and dim is the dimension number of each hawk.The mutation operator mathematical model of the proposed hybrid strategies can be expressed as(28)d=m∗UBj−LBj,j∈rand1,2,…,dim,Xt+1=Xj=Xj+d,where m is a random value in scope of (−0.5, 0.5).The specific pseudocode of the proposed HMOHHO algorithm is depicted in Algorithm1.Algorithm 1: Pseudocode of HMOHHO algorithm.
Inputs: Population size N, maximum iterations T and archive size SInitialize: Latin hypercube sampling initialize population Xii=1,2,…,NCalculate f1(·) and f2(·)Find the non-dominated solutionswhilet < Tifthe archive is fullCrowded area of the archive is deleted by roulette wheel methodAdd new solutions to the archiveend ifSelect Xrabbit(t) in the archiveforeach hawkifabs (E) < 1X(t) is updated by (14)elseX(t) is updated by (15)–(24)end ifXr1 is updated by (27)iff(Xr1) dominate f(X(t))X(t) = Xr1end ifXr2 is updated by (28)iff(Xr2) dominate f(X(t))X(t) = Xr2end ifend forCalculate f1(·) and f2(·)Find the non-dominated solutionst = t + 1end whilereturnarchive


## 4.2.1. MOHHO

Multiobjective Harris hawks optimization algorithm (MOHHO) based on HHO is expected to solve multiobjective optimization problem. The steps of MOHHO algorithm are as follows:Step 1: preset the population sizeN, maximum iterations T, and archive size S at the beginning as well as initializing the population randomly within the given ranges of variablesStep 2: calculate the objective valuesf1(·) and f2(·) for each hawk X(t)Step 3: select the nondominated Pareto-optimal solutions in the population at the current iteration; thus, the position of rabbitXrabbit(t) is acquired by leader selection mechanism, whose role is taking advantage of crowding distance to select a solution through roulette wheel method from a less populated area of the archiveStep 4: update each hawkX(t) according to equations (14)–(24)Step 5: calculate the new objective value for each hawk and find the nondominated solutions; thus better solutions will be recorded in the archiveStep 7: when the archive is full, the crowded area of the archive is deleted by the roulette wheel method for adding new solutions to the archiveStep 8: output the archive solutions


## 4.2.2. Reinforced MOHHO

Previous literatures [51–54] have shown that heuristic search algorithms generally suffer from the problem of trapping in local optimum easily, leading to the loss of solution diversity. Similarly, the original iterative process of MOHHO is not sufficient to maintain lateral diversity and achieve a Pareto front of good astringency and high diversity. In order to promote the searching capacity of MOHHO, the linear model of rabbit energy is replaced with a nonlinear one, as shown in(25)E=2E01−tT6.In addition to the improvement above, hybrid strategies are merged into MOHHO to jump out local optimum and search for more no-domain solutions. Firstly, population initialization is conducted through Latin hypercube sampling, which is the latest development in sampling technology [55]. The steps of Latin hypercube sampling initialization are as follows: (1) divide each dimension into N intervals that do not overlap each other; (2) randomly select a point within each interval for each dimension; (3) combine them into a vector. The mathematical model can be described by the following formula:(26)xd_min=xd0<xd1<xd2<⋯<xdj<⋯<xdN=xd_max,where xd donates the d-th dimension of each hawk vector and N donates population size.The modified differential evolution operator of the proposed hybrid strategies is shown in(27)Xt+1=Xrand1+F∗Xrand2−Xrand3,where Xrand(1), Xrand(2), and Xrand(3) are the positions of hawks chosen randomly in the current population, F is a 1 × dim vector whose value intervals for all dimensions are 0.3 and 0.8, and dim is the dimension number of each hawk.The mutation operator mathematical model of the proposed hybrid strategies can be expressed as(28)d=m∗UBj−LBj,j∈rand1,2,…,dim,Xt+1=Xj=Xj+d,where m is a random value in scope of (−0.5, 0.5).The specific pseudocode of the proposed HMOHHO algorithm is depicted in Algorithm1.Algorithm 1: Pseudocode of HMOHHO algorithm.
Inputs: Population size N, maximum iterations T and archive size SInitialize: Latin hypercube sampling initialize population Xii=1,2,…,NCalculate f1(·) and f2(·)Find the non-dominated solutionswhilet < Tifthe archive is fullCrowded area of the archive is deleted by roulette wheel methodAdd new solutions to the archiveend ifSelect Xrabbit(t) in the archiveforeach hawkifabs (E) < 1X(t) is updated by (14)elseX(t) is updated by (15)–(24)end ifXr1 is updated by (27)iff(Xr1) dominate f(X(t))X(t) = Xr1end ifXr2 is updated by (28)iff(Xr2) dominate f(X(t))X(t) = Xr2end ifend forCalculate f1(·) and f2(·)Find the non-dominated solutionst = t + 1end whilereturnarchive


## 4.3. Computational Complexity Analysis

The complexity analysis of each iteration in the proposed HMOHHO is as follows: assume that individual number isN and the dimension of each hawk is dim, and the number of individuals in the current archive set is assumed to be A. In HMOHHO, the rabbit position in archive set is selected by leader selection mechanism, during which process the computational complexity can be denoted by O(A). Then, computational complexity of location updating can be denoted by O(N∗ dim). Additionally, the computational complexity of the proposed hybrid strategies can be denoted by O(2 ∗N∗ dim) and the computational complexity of updating archive set can be denoted by O(A ∗N∗ dim). Next, the computational complexity of grouping and sorting archive sets is O(A2). On the whole, the computational complexity of each iteration in HMOHHO is O(A2) +O(A∗N∗ dim).


## 5. Algorithm Performance Verification

### 5.1. Verification Experiment Design

In order to verify the availability of the proposed HMOHHO algorithm, four multimode UF test functions and three ZDT test functions [56, 57] are employed to prove the performance. The test functions are shown in Table 2. Besides, NSGA-III, MOPSO, MOGWO, and MOHHO algorithms are introduced for comparison.Table 2
Test function problems.
NameMathematical formulationUF1f1=x1+2/J1∑j∈J1xj−sin6πx1+jπ/n2,f2=1−x+2/J2∑j∈J2xj−sin6πx1+jπ/n2J1=jjis odd and2≤j≤n,J2=jjis even and2≤j≤nUF2f1=x1+2/J1∑j∈J1yi2,f2=1−x+2/J2∑j∈J2yi2J1=jjis odd and2≤j≤n,J2=jjis even and2≤j≤nyi=xj−0.3x12cos24πx1+4jπ/n+0.6x1cos6πx1+jπ/nifj∈J1xj−0.3x12cos24πx1+4jπ/n+0.6x1cos6πx1+jπ/nifj∈J2UF4f1=x1+2/J1∑j∈J1hyj,f2=1−x2+2/J2∑j∈J2hyjJ1andJ2are identical to those of UF1,yj=xj−sin6πx1+jπ/n,j=2,3,…,n,ht=t/1+e2tUF7f1=x15+2/J1∑j∈J1yj2,f2=x15+2/J2∑j∈J2yj2J1andJ2are identical to those of UF1,yi=xj−sin6πx1+jπ/n,j=2,3,…,nZDT1f1X=x1f2X=g∗1−f1/ggX=1+9∗∑i=2nxi/n−1ZDT3f1X=x1f2X=g∗1−f1/g−f1/gsin10πf1gX=1+9∗∑i=2nxi/n−1ZDT6f1X=1−exp−4x1sin66πx1f2X=g∗1−f1/g2gX=1+9∗∑i=2nxi/n−10.25The UF test functions search space is [0, 1] × [−1, 1]dim−1 and the ZDT test functions search space is [0, 1]dim.


### 5.2. Parameter Setting

The principle and implementation of different algorithms are not uniform exactly. Consequently, relatively fair comparisons of different algorithms are achieved by setting the same maximum iterations, individual dimensions of population, population size, and archive size. The above four parameters of all algorithms for UF problems are set to 1000, 30, 100, and 100 severally. For problems ZDT1 and ZDT3, the parameters are set to 200, 30, 100, and 100 severally. For problem ZDT6, the parameters are set to 200, 10, 100, and 100, respectively. In MOGWO and MOPSO, the values of inflation rate, leader selection pressure, deletion selection pressure, and number of grids per dimension are both set to 0.1, 4, 2, and 10, respectively. The values of inertia, inertia weight damping rate, and mutation rate are set to 0.5, 0.99, and 0.1 in MOPSO, respectively. The values of crossover rate, number of neighbors, and mutation rate in NSGAIII are set to 0.5, 10, and 0.1, respectively.


### 5.3. Assessment Metrics

To evaluate the algorithm performance, there exist a number of assessment metrics, such as generational distance (GD), spacing (SP), inverted generational distance (IGD), maximum spread (MS), hypervolume (HV), and diversity metric (Δ) [58]. Among the previous metrics, IGD and HV which can measure convergence and diversity of solutions are chosen to evaluate the performance of each algorithm.HV means the volume of region in objective space enclosed by the nondominated solution set. The larger value of HV means the better overall performance of algorithm. The HV assessment metrics is defined as(29)HVp,zr=VOL∪X∈pf1X,z1r×⋯fmXzmr,where zr=z1r,z2r,…,zmr is the point in objective space dominated by any Pareto advantage and VOL(·) is the Lebesgue measure used to evaluate the volume.IGD means the average of distances from each reference point to the nearest solution. The smaller value of IGD means the better overall performance of algorithm. The IGD indicator is defined as(30)IGDP∗,P=∑X∈p∗dX,pp∗,where P is the solution set obtained by UF problem, P∗ is a set composed of Pareto-optimal reference points which are uniform distribution in Pareto front, and d(X, p) is the Euclidean distance between X and any point in P.


### 5.4. Performance Analysis

Each algorithm is run 30 times independently to eliminate contingency. The HV and IGD statistical results of different algorithms evaluated on each problem are shown in Table3. The Pareto-optimal solutions optimized by HMOHHO are shown in Figures 5 and 6.Table 3
HV and IGD statistical results of different algorithms.
AlgorithmHVIGDMinMaxMeanStdMinMaxMeanStdUF1 problem (2 objs)NSGAIII0.19400.42900.36420.61800.32190.73990.43560.1137MOPSO0.47780.58430.54650.02740.08180.16600.11150.0197MOGWO0.35800.59840.55550.04410.09150.25770.11910.0280MOHHO0.52630.57760.52630.03250.09470.18160.12980.0219HMOHHO0.61980.64400.63270.00580.05330.07350.06190.0050UF2 problem (2 objs)NSGAIII0.39400.53710.47100.03630.18970.50960.35990.0791MOPSO0.62610.66260.64850.00870.04790.12790.06260.0147MOGWO0.61700.65770.64250.01010.05050.08630.06520.0084MOHHO0.61820.63690.62590.00460.06960.08760.08110.0048HMOHHO0.64820.67210.65800.00520.04120.05980.05280.0041UF4 problem (2 objs)NSGAIII0.06490.30200.21900.05920.09690.59740.23590.1008MOPSO0.21110.25540.22630.00950.13230.17780.16170.0095MOGWO0.31100.35930.35600.00400.05800.07910.06230.0054MOHHO0.34400.35650.35140.00380.05950.07060.06320.0038HMOHHO0.34650.36230.35690.00480.05640.06820.06080.0042UF7 problem (2 objs)NSGAIII0.11700.35900.26760.58700.25840.66330.40560.0999MOPSO0.11010.49430.40880.10380.06180.67390.16600.1541MOGWO0.29600.50500.48070.03570.05850.33930.07910.0495MOHHO0.37650.48400.43660.02690.06380.14580.09940.0204HMOHHO0.50940.53590.52640.00560.03370.05050.03940.0038ZDT1 problem (2 objs)NSGAIII0.30700.50690.44600.04940.18930.47660.26580.0729MOPSO0.06280.66670.40220.16860.04270.72230.29070.1837MOGWO0.69800.71360.70780.00320.00840.02710.01470.0037MOHHO0.71050.71530.71290.00110.00880.01290.01080.0010HMOHHO0.71300.71800.71620.00100.00610.01010.00779.19E − 04ZDT3 problem (2 objs)NSGAIII0.40500.63390.53700.05910.20590.37660.26850.0468MOPSO0.01570.77560.30200.18440.13590.93770.53300.1860MOGWO0.56600.59900.58100.00630.00860.02490.01350.0037MOHHO0.57770.59060.58090.00220.00810.01300.01060.0013HMOHHO0.57970.58970.58210.00180.00700.00990.00859.73E − 05ZDT6 problem (2 objs)NSGAIII0.16700.38330.35100.05080.00300.07670.01230.0153MOPSO0.31350.38620.37710.02130.00260.02080.00570.0048MOGWO0.36500.38570.37980.00580.00330.10360.02460.0284MOHHO0.38160.38530.38370.00080.00450.02150.00650.0030HMOHHO0.38460.38660.38575.75E− 040.00290.00650.00428.27E − 04Figure 5
Pareto-optimal solutions obtained by HMOHHO for UF1, UF2, UF4, and UF7. (a) HMOHHO for UF1 problem, (b) HMOHHO for UF2 problem, (c) HMOHHO for UF4 problem, and (d) HMOHHO for UF7 problem.
(a)(b)(c)(d)Figure 6
Pareto-optimal solutions obtained by HMOHHO for ZDT1, ZDT3, and ZDT6. (a) HMOHHO for ZDT1 problem, (b) HMOHHO for ZDT3 problem, and (c) HMOHHO for ZDT6 problem.
(a)(b)(c)It can be seen from Table3 that the proposed HMOHHO algorithm achieves the best results on all the statistical indexes for UF1, UF7, and ZDT1, which means that the proposed HMOHHO algorithm provides superior diversity and astringency on UF1, UF7, and ZDT1, and the front coverage of Pareto-optimal solutions of HMOHHO is broader than that of all other algorithms on these test functions.As to the test functions UF2 and UF4, the proposed HMOHHO algorithm achieves the best HV and IGD results in terms of min, max, and mean values. Although MOHHO obtains the minimum standard deviation of HV for both UF2 and UF4 problems as well as the minimum standard deviation of IGD for UF4 problem, the standard deviation values of the proposed HMOHHO algorithm are close to the minimum standard deviations. On the whole, HMOHHO achieves the best convergence average with well stability. It also can be seen that the worst results belong to NSGAIII due to the fact that NSGAIII easily falls into local optimums.As to the test function ZDT3, the proposed HMOHHO algorithm achieves the best HV results in terms of min, mean, and std values as well as the best IGD in terms of all four values. Although MOPSO obtains the maximum HV value, its HV results in terms of mean and std values are not so satisfactory. With regard to the test function ZDT6, the proposed HMOHHO algorithm achieves the best IGD results in terms of max, mean, and std values as well as the best HV in terms of all four values. Although MOPSO obtains the minimum IGD value, its IGD results in terms of mean and std values are not so satisfactory.On the whole, based on the above contrastive analysis for Table3, the application for test functions shows that the performance of the proposed HMOHHO algorithm achieves better convergence and stability than other contrastive algorithms, indicating that the proposed HMOHHO is able to realize remarkable diversity and astringency ability in settling multiobjective problems.


## 5.1. Verification Experiment Design

In order to verify the availability of the proposed HMOHHO algorithm, four multimode UF test functions and three ZDT test functions [56, 57] are employed to prove the performance. The test functions are shown in Table 2. Besides, NSGA-III, MOPSO, MOGWO, and MOHHO algorithms are introduced for comparison.Table 2
Test function problems.
NameMathematical formulationUF1f1=x1+2/J1∑j∈J1xj−sin6πx1+jπ/n2,f2=1−x+2/J2∑j∈J2xj−sin6πx1+jπ/n2J1=jjis odd and2≤j≤n,J2=jjis even and2≤j≤nUF2f1=x1+2/J1∑j∈J1yi2,f2=1−x+2/J2∑j∈J2yi2J1=jjis odd and2≤j≤n,J2=jjis even and2≤j≤nyi=xj−0.3x12cos24πx1+4jπ/n+0.6x1cos6πx1+jπ/nifj∈J1xj−0.3x12cos24πx1+4jπ/n+0.6x1cos6πx1+jπ/nifj∈J2UF4f1=x1+2/J1∑j∈J1hyj,f2=1−x2+2/J2∑j∈J2hyjJ1andJ2are identical to those of UF1,yj=xj−sin6πx1+jπ/n,j=2,3,…,n,ht=t/1+e2tUF7f1=x15+2/J1∑j∈J1yj2,f2=x15+2/J2∑j∈J2yj2J1andJ2are identical to those of UF1,yi=xj−sin6πx1+jπ/n,j=2,3,…,nZDT1f1X=x1f2X=g∗1−f1/ggX=1+9∗∑i=2nxi/n−1ZDT3f1X=x1f2X=g∗1−f1/g−f1/gsin10πf1gX=1+9∗∑i=2nxi/n−1ZDT6f1X=1−exp−4x1sin66πx1f2X=g∗1−f1/g2gX=1+9∗∑i=2nxi/n−10.25The UF test functions search space is [0, 1] × [−1, 1]dim−1 and the ZDT test functions search space is [0, 1]dim.


## 5.2. Parameter Setting

The principle and implementation of different algorithms are not uniform exactly. Consequently, relatively fair comparisons of different algorithms are achieved by setting the same maximum iterations, individual dimensions of population, population size, and archive size. The above four parameters of all algorithms for UF problems are set to 1000, 30, 100, and 100 severally. For problems ZDT1 and ZDT3, the parameters are set to 200, 30, 100, and 100 severally. For problem ZDT6, the parameters are set to 200, 10, 100, and 100, respectively. In MOGWO and MOPSO, the values of inflation rate, leader selection pressure, deletion selection pressure, and number of grids per dimension are both set to 0.1, 4, 2, and 10, respectively. The values of inertia, inertia weight damping rate, and mutation rate are set to 0.5, 0.99, and 0.1 in MOPSO, respectively. The values of crossover rate, number of neighbors, and mutation rate in NSGAIII are set to 0.5, 10, and 0.1, respectively.


## 5.3. Assessment Metrics

To evaluate the algorithm performance, there exist a number of assessment metrics, such as generational distance (GD), spacing (SP), inverted generational distance (IGD), maximum spread (MS), hypervolume (HV), and diversity metric (Δ) [58]. Among the previous metrics, IGD and HV which can measure convergence and diversity of solutions are chosen to evaluate the performance of each algorithm.HV means the volume of region in objective space enclosed by the nondominated solution set. The larger value of HV means the better overall performance of algorithm. The HV assessment metrics is defined as(29)HVp,zr=VOL∪X∈pf1X,z1r×⋯fmXzmr,where zr=z1r,z2r,…,zmr is the point in objective space dominated by any Pareto advantage and VOL(·) is the Lebesgue measure used to evaluate the volume.IGD means the average of distances from each reference point to the nearest solution. The smaller value of IGD means the better overall performance of algorithm. The IGD indicator is defined as(30)IGDP∗,P=∑X∈p∗dX,pp∗,where P is the solution set obtained by UF problem, P∗ is a set composed of Pareto-optimal reference points which are uniform distribution in Pareto front, and d(X, p) is the Euclidean distance between X and any point in P.


## 5.4. Performance Analysis

Each algorithm is run 30 times independently to eliminate contingency. The HV and IGD statistical results of different algorithms evaluated on each problem are shown in Table3. The Pareto-optimal solutions optimized by HMOHHO are shown in Figures 5 and 6.Table 3
HV and IGD statistical results of different algorithms.
AlgorithmHVIGDMinMaxMeanStdMinMaxMeanStdUF1 problem (2 objs)NSGAIII0.19400.42900.36420.61800.32190.73990.43560.1137MOPSO0.47780.58430.54650.02740.08180.16600.11150.0197MOGWO0.35800.59840.55550.04410.09150.25770.11910.0280MOHHO0.52630.57760.52630.03250.09470.18160.12980.0219HMOHHO0.61980.64400.63270.00580.05330.07350.06190.0050UF2 problem (2 objs)NSGAIII0.39400.53710.47100.03630.18970.50960.35990.0791MOPSO0.62610.66260.64850.00870.04790.12790.06260.0147MOGWO0.61700.65770.64250.01010.05050.08630.06520.0084MOHHO0.61820.63690.62590.00460.06960.08760.08110.0048HMOHHO0.64820.67210.65800.00520.04120.05980.05280.0041UF4 problem (2 objs)NSGAIII0.06490.30200.21900.05920.09690.59740.23590.1008MOPSO0.21110.25540.22630.00950.13230.17780.16170.0095MOGWO0.31100.35930.35600.00400.05800.07910.06230.0054MOHHO0.34400.35650.35140.00380.05950.07060.06320.0038HMOHHO0.34650.36230.35690.00480.05640.06820.06080.0042UF7 problem (2 objs)NSGAIII0.11700.35900.26760.58700.25840.66330.40560.0999MOPSO0.11010.49430.40880.10380.06180.67390.16600.1541MOGWO0.29600.50500.48070.03570.05850.33930.07910.0495MOHHO0.37650.48400.43660.02690.06380.14580.09940.0204HMOHHO0.50940.53590.52640.00560.03370.05050.03940.0038ZDT1 problem (2 objs)NSGAIII0.30700.50690.44600.04940.18930.47660.26580.0729MOPSO0.06280.66670.40220.16860.04270.72230.29070.1837MOGWO0.69800.71360.70780.00320.00840.02710.01470.0037MOHHO0.71050.71530.71290.00110.00880.01290.01080.0010HMOHHO0.71300.71800.71620.00100.00610.01010.00779.19E − 04ZDT3 problem (2 objs)NSGAIII0.40500.63390.53700.05910.20590.37660.26850.0468MOPSO0.01570.77560.30200.18440.13590.93770.53300.1860MOGWO0.56600.59900.58100.00630.00860.02490.01350.0037MOHHO0.57770.59060.58090.00220.00810.01300.01060.0013HMOHHO0.57970.58970.58210.00180.00700.00990.00859.73E − 05ZDT6 problem (2 objs)NSGAIII0.16700.38330.35100.05080.00300.07670.01230.0153MOPSO0.31350.38620.37710.02130.00260.02080.00570.0048MOGWO0.36500.38570.37980.00580.00330.10360.02460.0284MOHHO0.38160.38530.38370.00080.00450.02150.00650.0030HMOHHO0.38460.38660.38575.75E− 040.00290.00650.00428.27E − 04Figure 5
Pareto-optimal solutions obtained by HMOHHO for UF1, UF2, UF4, and UF7. (a) HMOHHO for UF1 problem, (b) HMOHHO for UF2 problem, (c) HMOHHO for UF4 problem, and (d) HMOHHO for UF7 problem.
(a)(b)(c)(d)Figure 6
Pareto-optimal solutions obtained by HMOHHO for ZDT1, ZDT3, and ZDT6. (a) HMOHHO for ZDT1 problem, (b) HMOHHO for ZDT3 problem, and (c) HMOHHO for ZDT6 problem.
(a)(b)(c)It can be seen from Table3 that the proposed HMOHHO algorithm achieves the best results on all the statistical indexes for UF1, UF7, and ZDT1, which means that the proposed HMOHHO algorithm provides superior diversity and astringency on UF1, UF7, and ZDT1, and the front coverage of Pareto-optimal solutions of HMOHHO is broader than that of all other algorithms on these test functions.As to the test functions UF2 and UF4, the proposed HMOHHO algorithm achieves the best HV and IGD results in terms of min, max, and mean values. Although MOHHO obtains the minimum standard deviation of HV for both UF2 and UF4 problems as well as the minimum standard deviation of IGD for UF4 problem, the standard deviation values of the proposed HMOHHO algorithm are close to the minimum standard deviations. On the whole, HMOHHO achieves the best convergence average with well stability. It also can be seen that the worst results belong to NSGAIII due to the fact that NSGAIII easily falls into local optimums.As to the test function ZDT3, the proposed HMOHHO algorithm achieves the best HV results in terms of min, mean, and std values as well as the best IGD in terms of all four values. Although MOPSO obtains the maximum HV value, its HV results in terms of mean and std values are not so satisfactory. With regard to the test function ZDT6, the proposed HMOHHO algorithm achieves the best IGD results in terms of max, mean, and std values as well as the best HV in terms of all four values. Although MOPSO obtains the minimum IGD value, its IGD results in terms of mean and std values are not so satisfactory.On the whole, based on the above contrastive analysis for Table3, the application for test functions shows that the performance of the proposed HMOHHO algorithm achieves better convergence and stability than other contrastive algorithms, indicating that the proposed HMOHHO is able to realize remarkable diversity and astringency ability in settling multiobjective problems.


## 6. Case Experiment

### 6.1. Comparison of Control Effects between Two Controllers

In this experiment, FOPID and PID controllers of HTGS were applied to obtain dynamic performance under no-load and on-load disturbance conditions. Experiment was run over a limited time frame of 20 s. Considering two objective functionsf1(·) and f2(·), the Pareto fronts obtained by HMOHHO for HTGS multiobjective control problem based on PID and FOPID controllers are shown in Figure 7. As seen in Figure 7, the Pareto front based on FOPID controller is completely located in the below portion of the result based on PID controller, indicating that the values of f1(·) and f2(·) are smaller based on FOPID controller. Thus, better parameters are deduced based on FOPID controller to achieve better dynamic performance for HTGS.Figure 7
Pareto fronts optimized by HMOHHO for multiobjective control problems of HTGS based on PID and FOPID controllers.Seven control schemes selected from the Pareto-optimal solution set are shown in Table4, where the first six ones are from the Pareto front with FOPID controller and the last one is from the Pareto front with PID controller. The results obtained by control schemes 3 and 7 in Table 4 are shown in Figure 8. In the case of no-load disturbance condition, the system with FOPID controller has a less response rise time and low overshoot, while the system with PID controller takes a long time to restore. Under on-load disturbance condition, the system with FOPID controller can recover quickly.Table 4
Control schemes selected from the Pareto-optimal solution set.
Control schemeParameters of HGTSPerformance indexesKpKiKdλμITAE under no-loadITAE under on-load15.65405.41582.29570.99410.30830.29530.024725.72115.00531.54680.99160.28260.25620.026535.71434.87781.96150.99060.26790.25150.027045.89884.86681.54680.99160.27660.24230.028655.81634.76061.52470.98910.21340.23550.031866.11254.19031.15520.98940.29180.21530.033077.37913.55820.68840.88160.0628Figure 8
Comparison of PID and FOPID under two operating conditions. (a) No-load operating condition; (b) on-load operating condition.
(a)(b)


### 6.2. The Pareto-Optimal Sets with Different Algorithms and Analysis

In this experiment, NSGA-III, MOPSO, MOGWO, MOHHO, and the proposed HMOHHO algorithms were applied to optimize HTGS under no-load and on-load operating conditions. The parameters of HTGS model are listed in Table1. The simulation time is 20 s. The maximum iteration is 150 and each algorithm is run 10 times independently. The best results for each algorithm in the repeated experiments are shown in Figure 9, from which it can be seen that the solution set of the proposed HMOHHO achieves more diversity than the solutions of other algorithms.Figure 9
The Pareto fronts of different algorithms. (a) The Pareto fronts of HMOHHO and NSGAIII. (b) The Pareto fronts of HMOHHO and MOPSO. (c) The Pareto fronts of HMOHHO and MOGWO. (d) The Pareto fronts of HMOHHO and MOHHO.
(a)(b)(c)(d)Three control schemes from the optimal solution set in Pareto front as shown in Table4 are selected to conduct a comparative experiment of transient process response. The transient process responses with the three control schemes are shown in Figure 10 and the corresponding performance indexes are listed in Table 4. As exhibited in Figure 10, all of the three control schemes realize good control stability under no-load and on-load disturbance conditions. Besides, it can be concluded from Table 4 that the ITAE indexes under no-load and on-load are with a certain contradiction. Specifically, scheme 5 achieves the smallest overshoot among the three schemes under no-load disturbance condition, while its stability is the worst under on-load disturbance condition. On the contrary, scheme 3 is the most stable under on-load disturbance condition, while the overshoot is larger than the other two control schemes under no-load disturbance condition. Scheme 4 is a balanced one which is suitable for the transformation of operating conditions. In order to achieve better control quality of HTGS, the control schemes with good control effect should be chosen by decision makers in the final solution according to the special requirements of HTGS.Figure 10
Transient process responses with different control schemes. (a) No-load operating condition; (b) on-load operating condition.
(a)(b)


## 6.1. Comparison of Control Effects between Two Controllers

In this experiment, FOPID and PID controllers of HTGS were applied to obtain dynamic performance under no-load and on-load disturbance conditions. Experiment was run over a limited time frame of 20 s. Considering two objective functionsf1(·) and f2(·), the Pareto fronts obtained by HMOHHO for HTGS multiobjective control problem based on PID and FOPID controllers are shown in Figure 7. As seen in Figure 7, the Pareto front based on FOPID controller is completely located in the below portion of the result based on PID controller, indicating that the values of f1(·) and f2(·) are smaller based on FOPID controller. Thus, better parameters are deduced based on FOPID controller to achieve better dynamic performance for HTGS.Figure 7
Pareto fronts optimized by HMOHHO for multiobjective control problems of HTGS based on PID and FOPID controllers.Seven control schemes selected from the Pareto-optimal solution set are shown in Table4, where the first six ones are from the Pareto front with FOPID controller and the last one is from the Pareto front with PID controller. The results obtained by control schemes 3 and 7 in Table 4 are shown in Figure 8. In the case of no-load disturbance condition, the system with FOPID controller has a less response rise time and low overshoot, while the system with PID controller takes a long time to restore. Under on-load disturbance condition, the system with FOPID controller can recover quickly.Table 4
Control schemes selected from the Pareto-optimal solution set.
Control schemeParameters of HGTSPerformance indexesKpKiKdλμITAE under no-loadITAE under on-load15.65405.41582.29570.99410.30830.29530.024725.72115.00531.54680.99160.28260.25620.026535.71434.87781.96150.99060.26790.25150.027045.89884.86681.54680.99160.27660.24230.028655.81634.76061.52470.98910.21340.23550.031866.11254.19031.15520.98940.29180.21530.033077.37913.55820.68840.88160.0628Figure 8
Comparison of PID and FOPID under two operating conditions. (a) No-load operating condition; (b) on-load operating condition.
(a)(b)


## 6.2. The Pareto-Optimal Sets with Different Algorithms and Analysis

In this experiment, NSGA-III, MOPSO, MOGWO, MOHHO, and the proposed HMOHHO algorithms were applied to optimize HTGS under no-load and on-load operating conditions. The parameters of HTGS model are listed in Table1. The simulation time is 20 s. The maximum iteration is 150 and each algorithm is run 10 times independently. The best results for each algorithm in the repeated experiments are shown in Figure 9, from which it can be seen that the solution set of the proposed HMOHHO achieves more diversity than the solutions of other algorithms.Figure 9
The Pareto fronts of different algorithms. (a) The Pareto fronts of HMOHHO and NSGAIII. (b) The Pareto fronts of HMOHHO and MOPSO. (c) The Pareto fronts of HMOHHO and MOGWO. (d) The Pareto fronts of HMOHHO and MOHHO.
(a)(b)(c)(d)Three control schemes from the optimal solution set in Pareto front as shown in Table4 are selected to conduct a comparative experiment of transient process response. The transient process responses with the three control schemes are shown in Figure 10 and the corresponding performance indexes are listed in Table 4. As exhibited in Figure 10, all of the three control schemes realize good control stability under no-load and on-load disturbance conditions. Besides, it can be concluded from Table 4 that the ITAE indexes under no-load and on-load are with a certain contradiction. Specifically, scheme 5 achieves the smallest overshoot among the three schemes under no-load disturbance condition, while its stability is the worst under on-load disturbance condition. On the contrary, scheme 3 is the most stable under on-load disturbance condition, while the overshoot is larger than the other two control schemes under no-load disturbance condition. Scheme 4 is a balanced one which is suitable for the transformation of operating conditions. In order to achieve better control quality of HTGS, the control schemes with good control effect should be chosen by decision makers in the final solution according to the special requirements of HTGS.Figure 10
Transient process responses with different control schemes. (a) No-load operating condition; (b) on-load operating condition.
(a)(b)


## 7. Conclusions

To enhance the controlling parameter tuning applicability of HTGS under changeable operating conditions, the multiobjective optimization problem of FOPID controller is built by considering no-load disturbance and on-load disturbance operation conditions, where the ITAE performance indicators under both operation conditions are employed as the objective functions. Then, the newly proposed Harris hawks optimization is extended to the multiobjective version MOHHO for solving the optimization problem. Additionally, the global searching capability of MOHHO is enhanced greatly through coupling MOHHO with hybrid strategies (HMOHHO) which include Latin hypercube sampling initialization, modified differential evolution operator, and mutation operator as well as replacing the linear model of rabbit energy with a nonlinear one. Subsequently, the proposed HMOHHO algorithm is tested on several test functions and compared with NSGAIII, MOGWO, MOPSO, and MOHHO, which verifies the effectiveness of the proposed HMOHHO algorithm. The practical application as well as contrastive analysis shows that the constructed multiobjective problem of FOPID controller is suitable for HTGS under changeable operating conditions and the proposed HMOHHO algorithm possesses outstanding feasibility and superiority among all the algorithms.



---

*Source: 9274980-2020-07-18.xml*