# Multiobjective Optimal Control for Hydraulic Turbine Governing System Based on an Improved MOGWO Algorithm

Authors: Xin Xia; Jie Ji; Chao-shun Li; Xiaoming Xue; Xiaolu Wang; Chu Zhang
Journal: Complexity
 (2019)

---

## Abstract

Hydraulic turbine governing system (HTGS) is essential equipment which regulates frequency and power of the power grids. In previous studies, optimal control of HTGS is always aiming at one single operation condition. The variation of operation conditions of HTGS is seldom considered. In this paper, multiobjective optimal function is proposed for HTGS under multiple operation conditions. In order to optimize the solution to the multiobjective problems, a novel multiobjective grey wolf optimizer algorithm with searching factor (sMOGWO) is also proposed with two improvements: adding searching step to search more no-domain solutions nearby the wolves and adjusting control parameters to keep exploration ability in later period. At first, the searching ability of the sMOGWO has been verified on several UF test problems by statistical analysis. And then, the sMOGWO is applied to optimize the solutions of the multiobjective problems of HTGS, while different algorithms are employed for comparison. The experimental results indicate that the sMOGWO is more effective algorithm and improves the control quality of the HTGS under multiple operation conditions.


---

## Body

## 1. Introduction

With the increase of people’s consciousness of environmental protection, more and more renewable energy has been applied to replace the traditional energy, such as wind power [1, 2] and solar power [3]. As wind power and solar power are connected with the large-scale power grid, people want to maximally utilize renewable energy on the premise that power balance is maintained [4–6]. However, the output of wind power and solar power is not steady due to the intermittent and fluctuating nature. Therefore the power grid needs regulating equipment to keep power balance. As the hydroelectric unit is easier for changing the output power, it has been widely used as power and frequency regulating equipment in power grid [7, 8].Hydraulic turbine governing system (HTGS) which mostly contains PID controller is important automatic control equipment of hydroelectric unit [9]. The HTGS can adjust the output power of hydroelectric unit by changing the opening of guide. Therefore the control quality of HTGS influences the regulation quality of hydroelectric unit directly.Many research works have been done to improve the control quality; these works can be mainly divided into two categories: proposed new control models instead of PID controller and optimized control parameters of HTGS. In [10], a fractional order PID controller, whose order of derivative portion and integral portion is not integer, is proposed for HTGS. The fractional order PID controller provides more flexibility in achieving control objective. In [11], a fuzzy-PID controller is designed to improve control quality of HTGS. In [12], the sliding mode variable structure control strategy led to HTGS which has better robustness and adaptability than PID controller. In [13], a new adaptive inverse control method based on the learning characteristic of neural network was proposed for HTGS to improve the dynamic and stationary performance. All the above control models are useful, but they ignore that it is hard to change the control device to reach the control models for most operating hydroelectric units. Thus control parameter optimization seems more suitable to improve control quality for actual situation. To optimize PID controller, some popular optimization algorithms, including genetic algorithm (GA) [14], particle swarm optimization (PSO) [15], gravitational search algorithm (GSA) [16, 17], and ant lion optimization (ALO) [18], have been successfully applied in parameter optimization of HTGS. All these algorithms optimize the PID parameters by optimizing one objective function which means the optimal parameters are under one single operating condition; few of the researches have taken multiple operation conditions into consideration. The optimal results may not be suitable for variation of operation conditions of HTGS. Therefore, the study on parameter optimal control of HTGS under multiple operation conditions is meaningful and challengeable.Parameter optimal control of HTGS under multiple operation conditions actually refers to the multiobjective optimization problems. Some multiobjective optimization algorithms have been proposed and successfully applied in various applications in decades [19–25], such as nondominated sorting genetic algorithm II (NSGA-II) [19], multiobjective evolutionary algorithm based on decomposition (MOEA/D) [21], the strength Pareto evolutionary algorithm (PESA-II) [22], multiobjective particle swarm optimization (MOPSO) [23], and the multiobjective grey wolf optimizer (MOGWO) [25]. Among these multiobjective optimization algorithms, MOGWO algorithm has proved to be effective in multiobjective optimization problems than other algorithms [25]. However, as the guidance of lead wolves is much greater than the random factor, this algorithm is easy to fall into local optimum and has poor stability. Therefore it is necessary to enhance the searching ability of MOGWO.The main contribution of this paper is the design and optimization of HTGS under multiple operation conditions by using an improved MOGWO algorithm. Firstly, the problem in optimal control of HTGS is explained and multiobjective optimal function is proposed. Secondly, a novel MOGWO algorithm based on searching factor (sMOGWO) is proposed to optimize the multiobjective problem. The sMOGWO is expected to improve the searching ability of MOGWO from two aspects: adding searching step to search more no-domain solutions nearby the wolves and adjusting control parameters to keep exploration ability in later period.The remaining part of this paper is organized as follows: the HTGS model and its control problem are discussed in Section2. In Section 3, the weakness of MOGWO algorithm is analyzed and an improved MOGWO algorithm is proposed. In Section 4, simulation verification is present to demonstrate the advantage of the proposed sMOGWO algorithm compared with MOGWO, MOPSO, and other popular multiobjective evolutionary algorithms. Section 5 illustrates the case study and results along with a few discussions. Finally, conclusions from this research as well as the advantages and limitations of the proposed algorithm are discussed in Section 6.


## 2. HTGS Model and Its Control Problems

### 2.1. The HTGS Model

The HTGS is essential equipment of the hydropower station. The HTGS contains PID controller, electrohydraulic servo system, hydroturbine, and hydrogenerator.The PID controller has been widely used in HTGS. In this paper, the structure of the PID controller model and electrohydraulic servo system is set as in Figure1. Many nonlinear factors have been considered in the electrohydraulic servo system, such as the dead zone and the relay device limiting, wherexr is the frequency giving;xu is the disturbance frequency;x is the frequency of the hydrogenerator unit;kp,kd,ki are the proportionality coefficient, differential coefficient, and integration coefficient, respectively; T1v is the differential filtering time constant;bp is the permanent difference coefficient;Ty is the response time constant of main control valve; andy represents the gate opening.Figure 1
PID control device and electric-hydraulic servo system.The hydroturbine is a complex system with strong nonlinear and time-varying character. Generally, the hydroturbine model can be descripted as follows:(1)mt=mty,h,ωq=qy,h,ωwheremt is the hydraulic machinery power,q is discharge of hydraulic turbine,y represents the gate opening, ω represents speed of hydraulic turbine, andh represents the water head.In this paper, the little fluctuation in transient of hydraulic power station is researched; the model of hydraulic turbine can be described as liner model, as follows:(2)mt=exω+eyy+ehhq=eqxω+eqyy+eqhhwhereeh is the transfer coefficient of turbine torque on the water head,ey is the transfer coefficient of turbine torque on the guide leaf opening,ex is the transfer coefficient of turbine torque on the speed, eqx is the transfer coefficient of turbine flow on the speed, eqy the transfer coefficient of turbine flow on guide leaf opening, and eqh is the transfer coefficient of turbine flow on the water head.In water diversion system, the water hammer and pipe wall can be considered as rigid under small fluctuations. The transfer function of rigid water diversion system can be expressed as(3)Grs=-Twswhere Tw is the water inertia time constant.In this paper, the core problem is the dynamic response of the HTGS. Hence, the dynamic speed of the generator is considered. The model of the generator can be described as the following transfer function:(4)xmt=1Tas+egwherex represents the frequency of the hydrogenerator unit,Ta represents the inertia time constant, and eg is the adaptive control coefficient.As in the above description of the mathematical model, the hydrogenerator unit model can be described as the block diagram in Figure2.Figure 2
Hydroturbine unit model with rigid water hammer.


### 2.2. Optimal Control of HTGS and Problems Analysis

The performance indexes of the optimal control of the HTGS are the most direct measure to evaluate the control quality. Therefore, the optimal control of the HTGS is the optimal control of performance indexes. The objective functions of optimal control are formed by a certain performance index or several performance indexes in most researches. The performance indexes in HTGS can be mainly divided into the following two categories [26]:(i) Performance indexes for the transition process mainly include adjusting time, rising time, steady state error, and overshoot.(ii) Performance indexes of error functional integral mainly include integral time absolute error (ITAE), integral absolute error (IAE), integral square time absolute error (ISTAE), integral square time square error (ISTSE), integral time square error (ITSE), and integral square error (ISE).The performance indexes of error functional integral are comprehensive indexes which make optimal control easier to meet the control requirements. The ITAE is one of the most widely used which has the characteristics of steady adjustment and small overshoot.Many researches have been done for optimal control of HTGS by optimizing ITAE. However most of the researches are for one single operation condition. As the hydrogenerator usually plays the roles of generating electricity, peak modulation, frequency modulation, and voltage modulation, there is a variety of operation conditions and frequent changes. The optimal control parameters may not be suitable for all the operation conditions.Here, we will give an example to explain this problem. Parameters of two classical operation conditions of the above HGTS model are given as follows:No-load:Ty1=0.02s,Ty=0.1s, Tw=0.66s, T1v=0.28,Ta=8.51s,bp=0.01, eg=1,eh=0.508,ey=0.903,ex=-0.242, eqx=0.634, eqy=-0.396, eqh=0.261.On-load:Ty1=0.02s,Ty=0.1s, Tw=0.66s, T1v=0.28, Ta=8.51s,bp=0, eg=1,eh=1.34,ey=0.926,ex=-1.184, eqx=0.371, eqy=0.974, eqh=0.308.5% step disturbance is set for HGTS. The optimal control parameters of HTGS under no-load condition are acquired by optimizing ITAE. The parameters areKp=14.665,Ki=3.026, andKd=0. The control transient process is shown in Figure 3. The control effect is well. However, if the parameters are used for the on-load condition, it is an ineffective control strategy.Figure 3
Transient process by different control strategies.The result indicates that the optimal control parameters of one single operation condition cannot be suitable for all the operation conditions. Accordingly, the optimal control of HTGS is a multiobjective optimization problem in fact. In this paper, we have proposed a novel control parameters optimization objective function. Different operation conditions are considered for optimal control. Because the no-load and on-load are the extreme operating conditions, the ITAE of both conditions are the objective functions. The objective function can be described as follows:(5)minJ1=ITAEno-load=f1Kp,Ki,Kd,yJ2=ITAEon-load=f2Kp,Ki,Kd,y(6)subjectto0≤Kp≤150≤Ki≤150≤Kd≤50≤y≤1wheref1 andf2 represent the relationship between the parameters and the ITAE. Equation (6) represents the range of parameters.


## 2.1. The HTGS Model

The HTGS is essential equipment of the hydropower station. The HTGS contains PID controller, electrohydraulic servo system, hydroturbine, and hydrogenerator.The PID controller has been widely used in HTGS. In this paper, the structure of the PID controller model and electrohydraulic servo system is set as in Figure1. Many nonlinear factors have been considered in the electrohydraulic servo system, such as the dead zone and the relay device limiting, wherexr is the frequency giving;xu is the disturbance frequency;x is the frequency of the hydrogenerator unit;kp,kd,ki are the proportionality coefficient, differential coefficient, and integration coefficient, respectively; T1v is the differential filtering time constant;bp is the permanent difference coefficient;Ty is the response time constant of main control valve; andy represents the gate opening.Figure 1
PID control device and electric-hydraulic servo system.The hydroturbine is a complex system with strong nonlinear and time-varying character. Generally, the hydroturbine model can be descripted as follows:(1)mt=mty,h,ωq=qy,h,ωwheremt is the hydraulic machinery power,q is discharge of hydraulic turbine,y represents the gate opening, ω represents speed of hydraulic turbine, andh represents the water head.In this paper, the little fluctuation in transient of hydraulic power station is researched; the model of hydraulic turbine can be described as liner model, as follows:(2)mt=exω+eyy+ehhq=eqxω+eqyy+eqhhwhereeh is the transfer coefficient of turbine torque on the water head,ey is the transfer coefficient of turbine torque on the guide leaf opening,ex is the transfer coefficient of turbine torque on the speed, eqx is the transfer coefficient of turbine flow on the speed, eqy the transfer coefficient of turbine flow on guide leaf opening, and eqh is the transfer coefficient of turbine flow on the water head.In water diversion system, the water hammer and pipe wall can be considered as rigid under small fluctuations. The transfer function of rigid water diversion system can be expressed as(3)Grs=-Twswhere Tw is the water inertia time constant.In this paper, the core problem is the dynamic response of the HTGS. Hence, the dynamic speed of the generator is considered. The model of the generator can be described as the following transfer function:(4)xmt=1Tas+egwherex represents the frequency of the hydrogenerator unit,Ta represents the inertia time constant, and eg is the adaptive control coefficient.As in the above description of the mathematical model, the hydrogenerator unit model can be described as the block diagram in Figure2.Figure 2
Hydroturbine unit model with rigid water hammer.


## 2.2. Optimal Control of HTGS and Problems Analysis

The performance indexes of the optimal control of the HTGS are the most direct measure to evaluate the control quality. Therefore, the optimal control of the HTGS is the optimal control of performance indexes. The objective functions of optimal control are formed by a certain performance index or several performance indexes in most researches. The performance indexes in HTGS can be mainly divided into the following two categories [26]:(i) Performance indexes for the transition process mainly include adjusting time, rising time, steady state error, and overshoot.(ii) Performance indexes of error functional integral mainly include integral time absolute error (ITAE), integral absolute error (IAE), integral square time absolute error (ISTAE), integral square time square error (ISTSE), integral time square error (ITSE), and integral square error (ISE).The performance indexes of error functional integral are comprehensive indexes which make optimal control easier to meet the control requirements. The ITAE is one of the most widely used which has the characteristics of steady adjustment and small overshoot.Many researches have been done for optimal control of HTGS by optimizing ITAE. However most of the researches are for one single operation condition. As the hydrogenerator usually plays the roles of generating electricity, peak modulation, frequency modulation, and voltage modulation, there is a variety of operation conditions and frequent changes. The optimal control parameters may not be suitable for all the operation conditions.Here, we will give an example to explain this problem. Parameters of two classical operation conditions of the above HGTS model are given as follows:No-load:Ty1=0.02s,Ty=0.1s, Tw=0.66s, T1v=0.28,Ta=8.51s,bp=0.01, eg=1,eh=0.508,ey=0.903,ex=-0.242, eqx=0.634, eqy=-0.396, eqh=0.261.On-load:Ty1=0.02s,Ty=0.1s, Tw=0.66s, T1v=0.28, Ta=8.51s,bp=0, eg=1,eh=1.34,ey=0.926,ex=-1.184, eqx=0.371, eqy=0.974, eqh=0.308.5% step disturbance is set for HGTS. The optimal control parameters of HTGS under no-load condition are acquired by optimizing ITAE. The parameters areKp=14.665,Ki=3.026, andKd=0. The control transient process is shown in Figure 3. The control effect is well. However, if the parameters are used for the on-load condition, it is an ineffective control strategy.Figure 3
Transient process by different control strategies.The result indicates that the optimal control parameters of one single operation condition cannot be suitable for all the operation conditions. Accordingly, the optimal control of HTGS is a multiobjective optimization problem in fact. In this paper, we have proposed a novel control parameters optimization objective function. Different operation conditions are considered for optimal control. Because the no-load and on-load are the extreme operating conditions, the ITAE of both conditions are the objective functions. The objective function can be described as follows:(5)minJ1=ITAEno-load=f1Kp,Ki,Kd,yJ2=ITAEon-load=f2Kp,Ki,Kd,y(6)subjectto0≤Kp≤150≤Ki≤150≤Kd≤50≤y≤1wheref1 andf2 represent the relationship between the parameters and the ITAE. Equation (6) represents the range of parameters.


## 3. An Improved MOGWO Algorithm

MOGWO algorithm as a new intelligent optimization algorithm has a better convergence rate than other intelligent optimization algorithms. Therefore, it has received great attention and wide application since it has been proposed.


### 3.1. Grey Wolf Optimizer and MOGWO Algorithm

The grey wolf optimizer (GWO) algorithm proposed by Mirjalili et al. in 2014 is a new intelligent optimization algorithm mimicking the hierarchies and hunting strategies of wolves [27]. In GWO, each grey wolf is treated as a potential solution. The grey wolves at the best solution, the second best solution, and the third best solution are treated as α, β, and δ wolves. The rest of the grey wolves are treated as ω wolves. The wolves are approaching the food position (global optimal solution) according to the guidance of α, β, and δ wolf. The equations that simulate the guidance of hunting are as follows:(7)D=C·Xpt-Xt(8)Xt+1=Xpt-A·DwhereX represents the position vector of a grey wolf,Xp represents the position vector of the prey, andt represents the current iteration.AandC are coefficient vectors subject to the following equations:(9)A=2·a·r1-a(10)C=2·r2wherer1 andr2 are random vectors from 0 to 1. a is control parameter from 2 to 0 and it will linearly decrease as the number of iterations increases.In 2016, Mirjalili proposed a multiobjective GWO algorithm named ‘MOGWO algorithm’ for multiobjective problems. There are two major changes in MOGWO: using the external population Archive to store the current nondominated solutions and proposing a selection strategy for multiobjective optimization.


### 3.2. sMOGWO Algorithm

Although the MOGWO algorithm has a better convergence rate, it is easy to fall into local optimum and has poor stability. The main reasons are summarized as follows:(i) The algorithm has great randomness only when initializing the position of wolves. Even though there is random factor in algorithm, when the position of wolves is updated, the effect of the guidance of lead wolves is much greater than the random factor. Thus the algorithm is highly dependent on the initial value, and self-regulation ability is weak.(ii) The MOGWO algorithms selectα, β, andδ wolves from Archive. However, if Archive set falls into local optimum, the algorithms can hardly skip the local optimum. Therefore, the exploration ability of the algorithm needs to be improved at later period of the optimization.As the guidance of lead wolves has too much influence, the grey wolves always blindly follow the lead wolves and the nondominated solutions around the lead wolves. The grey wolves always ignore the nondominated solutions beside them. Actually, the grey wolves often pass nearby other nondominant solutions when following the lead wolves. If the grey wolves have the ability of independent searching, the global optimization ability of the algorithm will be greatly improved.Consequently, in this paper, we proposed to add the searching factor in MOGWO named ‘sMOGWO algorithm’. After the position updating of wolves, the wolves will have the searching step. Each wolf will search the nearby position randomly; if the nearby position is better than the current position, the wolf will move to the new position. The step can be described as the following mathematical expressions:(11)Δx1,x2,⋯xdim=xk=r·ubk-lbk,k=rand1,2,⋯dimxi=0,i=1,2,⋯dimandi≠k(12)Xsearching=X+Δ(13)Xnew=Xsearching,iffXsearching≻fXXnew=X,otherwisewhere Δ is the searching step of the wolf,r is a random number between -0.5 and 0.5,ubk andlbk are the upper and lower boundary ofxk, dim is the dimension of wolf,X is the position of the wolf,Xnew is the new position after searching, and ≻ represents the dominance relationship.In the original algorithm, the control parametera linearly decreases as the number of iterations increases to keep high exploration ability at early period and high exploitation ability at later period. However, there is higher demands for exploration ability in multiobjective optimization. Therefore, in this paper, the control parametera is also modified. The adjustment strategy of control parametera must be subject to the following:(14)a=2-2·fx,x=itMaxit∈0.1f0=0,f1=1whereMaxit is the maximum number of iterations andf(x) is a monotone increasing function, in original algorithmf(x)=x.The higher the value ofa, the stronger the exploration ability of the algorithm, so thef(x) should be a concave function among [0,1]. Here, we proposed to use the following equation for control parametera.(15)a=2-2·itMaxit2The flow chart of the sMOGWO algorithm is shown in Figure 4.Figure 4
The flow chart of the sMOGWO algorithm.The steps of the sMOGWO algorithm are as follows.Step 1.
Archive size isNA, population size of grey wolves isN, number of grids per dimension isdim, inflation rate isa, and maximum number of iterations isMaxit. Leader selection pressure and deletion selection pressure are also set. The grey wolves and running parameters are initialized.Step 2.
The objective function values of each individual wolf are calculated. The Archive set is established, and the Archive population is grouped into iterative process.Step 3.
TheA→, C→, anda are calculated according to (9), (10), and (15); then the lead wolves are selected.Step 4.
The position of wolves is updated according to (7) and (8).Step 5.
The wolves get into the searching step according to (11), (12), and (13)Step 6.
The objective function values of each individual in wolves are calculated. The nondominant solutions of wolves are compared to the individuals in the Archive population one by one, and the Archive set is updated.Step 7.
The updated Archive set is regrouped and the number of individuals in the Archive population is checked. If the number of individuals exceeds the maximum of population, the extra solutions will be deleted.Step 8.
The algorithm ends when the maximum number of iterations is reached. All individuals in the Archive population are the optimization results of the algorithm. Otherwise, the algorithm returns back to Step3.


### 3.3. Computational Complexity Analysis

The computational complexity of MOGWO and sMOGWO is analyzed in this part.The computational complexity of each iteration (Step3 to Step 7) of the iterative process is analyzed one by one. The number of individuals in the current Archive set isA, the number of individuals in the grey wolf population isN, and the dimension of each individual isdim.In Step3, lead wolves are selected by calculating the roulette probability for each individual in the Archive set. The computational complexity can be expressed as O(A).In Step4, the computational complexity of position updating can be expressed as O(N∗dim).In Step5, the objective function values need be calculated by searching step, and the computational complexity can be expressed as O(N∗dim).In Step6, the nondominant solutions of wolves are compared to the individuals in the Archive set one by one to update the Archive set. The computational complexity is O(A∗N∗dim).In Step7, the Archive sets are grouped and sorted, and the computational complexity is O(A2).The computational complexity of sMOGWO in each iteration is O(A2)+ O(A∗N∗dim). The MOGWO includes all the above steps expect Step 5, and the computational complexity of MOGWO is also O(A2)+ O(A∗N∗dim). The proposed modified method has the same computational complexity as the original method. In addition, although an additional Step 5 is added to the sMOGWO algorithm, it can effectively reduce Archive set with a large number of similar nondominant solutions, which has great influence on the computational complexity. In fact, the proposed method can reduce the optimization time.


## 3.1. Grey Wolf Optimizer and MOGWO Algorithm

The grey wolf optimizer (GWO) algorithm proposed by Mirjalili et al. in 2014 is a new intelligent optimization algorithm mimicking the hierarchies and hunting strategies of wolves [27]. In GWO, each grey wolf is treated as a potential solution. The grey wolves at the best solution, the second best solution, and the third best solution are treated as α, β, and δ wolves. The rest of the grey wolves are treated as ω wolves. The wolves are approaching the food position (global optimal solution) according to the guidance of α, β, and δ wolf. The equations that simulate the guidance of hunting are as follows:(7)D=C·Xpt-Xt(8)Xt+1=Xpt-A·DwhereX represents the position vector of a grey wolf,Xp represents the position vector of the prey, andt represents the current iteration.AandC are coefficient vectors subject to the following equations:(9)A=2·a·r1-a(10)C=2·r2wherer1 andr2 are random vectors from 0 to 1. a is control parameter from 2 to 0 and it will linearly decrease as the number of iterations increases.In 2016, Mirjalili proposed a multiobjective GWO algorithm named ‘MOGWO algorithm’ for multiobjective problems. There are two major changes in MOGWO: using the external population Archive to store the current nondominated solutions and proposing a selection strategy for multiobjective optimization.


## 3.2. sMOGWO Algorithm

Although the MOGWO algorithm has a better convergence rate, it is easy to fall into local optimum and has poor stability. The main reasons are summarized as follows:(i) The algorithm has great randomness only when initializing the position of wolves. Even though there is random factor in algorithm, when the position of wolves is updated, the effect of the guidance of lead wolves is much greater than the random factor. Thus the algorithm is highly dependent on the initial value, and self-regulation ability is weak.(ii) The MOGWO algorithms selectα, β, andδ wolves from Archive. However, if Archive set falls into local optimum, the algorithms can hardly skip the local optimum. Therefore, the exploration ability of the algorithm needs to be improved at later period of the optimization.As the guidance of lead wolves has too much influence, the grey wolves always blindly follow the lead wolves and the nondominated solutions around the lead wolves. The grey wolves always ignore the nondominated solutions beside them. Actually, the grey wolves often pass nearby other nondominant solutions when following the lead wolves. If the grey wolves have the ability of independent searching, the global optimization ability of the algorithm will be greatly improved.Consequently, in this paper, we proposed to add the searching factor in MOGWO named ‘sMOGWO algorithm’. After the position updating of wolves, the wolves will have the searching step. Each wolf will search the nearby position randomly; if the nearby position is better than the current position, the wolf will move to the new position. The step can be described as the following mathematical expressions:(11)Δx1,x2,⋯xdim=xk=r·ubk-lbk,k=rand1,2,⋯dimxi=0,i=1,2,⋯dimandi≠k(12)Xsearching=X+Δ(13)Xnew=Xsearching,iffXsearching≻fXXnew=X,otherwisewhere Δ is the searching step of the wolf,r is a random number between -0.5 and 0.5,ubk andlbk are the upper and lower boundary ofxk, dim is the dimension of wolf,X is the position of the wolf,Xnew is the new position after searching, and ≻ represents the dominance relationship.In the original algorithm, the control parametera linearly decreases as the number of iterations increases to keep high exploration ability at early period and high exploitation ability at later period. However, there is higher demands for exploration ability in multiobjective optimization. Therefore, in this paper, the control parametera is also modified. The adjustment strategy of control parametera must be subject to the following:(14)a=2-2·fx,x=itMaxit∈0.1f0=0,f1=1whereMaxit is the maximum number of iterations andf(x) is a monotone increasing function, in original algorithmf(x)=x.The higher the value ofa, the stronger the exploration ability of the algorithm, so thef(x) should be a concave function among [0,1]. Here, we proposed to use the following equation for control parametera.(15)a=2-2·itMaxit2The flow chart of the sMOGWO algorithm is shown in Figure 4.Figure 4
The flow chart of the sMOGWO algorithm.The steps of the sMOGWO algorithm are as follows.Step 1.
Archive size isNA, population size of grey wolves isN, number of grids per dimension isdim, inflation rate isa, and maximum number of iterations isMaxit. Leader selection pressure and deletion selection pressure are also set. The grey wolves and running parameters are initialized.Step 2.
The objective function values of each individual wolf are calculated. The Archive set is established, and the Archive population is grouped into iterative process.Step 3.
TheA→, C→, anda are calculated according to (9), (10), and (15); then the lead wolves are selected.Step 4.
The position of wolves is updated according to (7) and (8).Step 5.
The wolves get into the searching step according to (11), (12), and (13)Step 6.
The objective function values of each individual in wolves are calculated. The nondominant solutions of wolves are compared to the individuals in the Archive population one by one, and the Archive set is updated.Step 7.
The updated Archive set is regrouped and the number of individuals in the Archive population is checked. If the number of individuals exceeds the maximum of population, the extra solutions will be deleted.Step 8.
The algorithm ends when the maximum number of iterations is reached. All individuals in the Archive population are the optimization results of the algorithm. Otherwise, the algorithm returns back to Step3.


## 3.3. Computational Complexity Analysis

The computational complexity of MOGWO and sMOGWO is analyzed in this part.The computational complexity of each iteration (Step3 to Step 7) of the iterative process is analyzed one by one. The number of individuals in the current Archive set isA, the number of individuals in the grey wolf population isN, and the dimension of each individual isdim.In Step3, lead wolves are selected by calculating the roulette probability for each individual in the Archive set. The computational complexity can be expressed as O(A).In Step4, the computational complexity of position updating can be expressed as O(N∗dim).In Step5, the objective function values need be calculated by searching step, and the computational complexity can be expressed as O(N∗dim).In Step6, the nondominant solutions of wolves are compared to the individuals in the Archive set one by one to update the Archive set. The computational complexity is O(A∗N∗dim).In Step7, the Archive sets are grouped and sorted, and the computational complexity is O(A2).The computational complexity of sMOGWO in each iteration is O(A2)+ O(A∗N∗dim). The MOGWO includes all the above steps expect Step 5, and the computational complexity of MOGWO is also O(A2)+ O(A∗N∗dim). The proposed modified method has the same computational complexity as the original method. In addition, although an additional Step 5 is added to the sMOGWO algorithm, it can effectively reduce Archive set with a large number of similar nondominant solutions, which has great influence on the computational complexity. In fact, the proposed method can reduce the optimization time.


## 4. Simulation Verification

In this section, simulation verification is present to demonstrate the advantage of the proposed sMOGWO algorithm. The proposed algorithm is compared to the MOGWO algorithm [25] and some other popular multiobjective optimal algorithms [28]: MOPSO, NSGA-III, MOEA/D, and SPEA2.


### 4.1. Test Problems

UF series test problems in CEC 2009 are multimode test functions [29], and there are a lot of locally optimal solutions, so it is suitable for test of multiobjective optimal algorithms. Three typical test problems are chosen as follows: UF3 whose Pareto boundary is concave shape, UF4 whose Pareto boundary is convex shape, UF7 whose Pareto boundary is line. The three UF problems are as follows:UF3 problem,n=30, the search space is 0,1n:(16)minf1=x1+2J14∑j∈J1yj2-2∏j∈J1cos⁡20yiπj+2f2=1-x1+2J14∑j∈J1yj2-2∏j∈J1cos⁡20yiπj+2subjecttoyi=xj-x10.51+3j-2/n-2,j=2,…,n,J1=j∣jisoddand2≤j≤nJ2=j∣jisevenand2≤j≤nUF4 problem,n=30, the search space is [0,1]×-2,2n-1:(17)minf1=x1+2J1∑j∈J1hyif2=1-x12+2J2∑j∈J1hyisubjecttoyi=xj-sin⁡6πx1+jπn,j=2,…,n,ht=t1+e2tJ1=j∣jisoddand2≤j≤nJ2=j∣jisevenand2≤j≤nUF7 problem,n=30, the search space is [0,1]×-1,1n-1:(18)minf1=x15+2J1∑j∈J1yi2f2=1-x15+2J2∑j∈J1yi2subjecttoyi=xj-sin⁡6πx1+jπn,j=2,…,n,J1=j∣jisoddand2≤j≤nJ2=j∣jisevenand2≤j≤n.


### 4.2. Parameters Setting

The principles and the operation modes of the algorithms are different. In order to make the algorithm contrast, the same maximum number of iterations, population size, and Archive size are set up. The parameters of different multiobjective optimization algorithms are shown in Table1.Table 1
The parameters of different multiobjective optimization algorithms.
Description sMOGWO MOGWO MOPSO NSGA-III MOEA/D SPEA2 Maximum number of iterations 1000 1000 1000 1000 1000 1000  Population size nPop=100 nPop=100 nPop=100 nPop=100 nPop=100 nPop=100  Archive size 100 100 100 100 100 100  Inflation rate α=0.1 α=0.1 α=0.1 / / /  Number of grids per dimension nGrid=10 nGrid=10 nGrid=10 / / /  Leader selection pressure β=4 β=4 β=4 / / /  Deletion selection pressure γ=2 γ=2 γ=2 / / /  Inertia weight / / ω=0.5 / / /  Inertia weight damping rate / / ω d a m p=0.99 / / /  Mutation rate / / mu=0.1 mu=0.1 mu=0.1 mu=0.1  Crossover rate / / / 0.5 0.5 0.5  Number of Neighbors / / / / 10 /


### 4.3. Performance Evaluation Indexes

The Generational Distance (GD) and the Spacing (SP) [30] have been used to evaluation the performance of each algorithm. (19)GD=1N∑i=1NDi2whereDi is the Euclidean distance between theith nondominant solution in Pareto solution set and the closest nondominant solution on the real Pareto front.N is the number of Pareto optimal solutions. The smaller value of GD indicates that the Pareto optimal solution is closer to the real Pareto front.(20)SP=df+dl+∑i=1N-1di-d-df+dl+N-1·d-wheredi is the Euclidean distance between the adjacent points of Pareto optimal solution.df anddl are the Euclidean distances between the endpoint of Pareto optimal solution and real Pareto front.N is the number of Pareto optimal solutions. The smaller value of SP indicates that Pareto optimal solution distribution is more homogenized.


### 4.4. Result Analysis

In order to eliminate the contingency, each algorithm is run 30 times independently. The statistical results of the evaluation indexes are shown in Tables2 and 3. In order to reflect the advantages and disadvantages of the optimization results more directly, the boxplots of GD value and SP value of each algorithm are given as in Figures 5 and 6.Table 2
Statistical results of GD for different algorithms.
GD sMOGWO MOGWO MOPSO NSGA-III MOEA/D SPEA2 UF3 Mean 0.0210 0.0542 0.0853 0.0268 0.0561 0.0254 Med 0.0211 0.0540 0.0737 0.0255 0.0571 0.0252 Max 0.0253 0.0900 0.2322 0.0716 0.0829 0.0342 Min 0.0180 0.0315 0.0425 0.0166 0.0299 0.0116  UF4 Mean 0.0239 0.0331 0.0361 0.0516 0.0409 0.0485 Med 0.0241 0.0330 0.0363 0.0521 0.0439 0.0475 Max 0.0257 0.0394 0.0406 0.0622 0.0526 0.0618 Min 0.0234 0.0259 0.0324 0.0415 0.0320 0.0340  UF7 Mean 0.0078 0.0267 0.0468 0.0319 0.0802 0.0313 Med 0.0077 0.0278 0.0424 0.0316 0.0584 0.0301 Max 0.0091 0.1012 0.0881 0.0726 0.1751 0.0608 Min 0.0072 0.0148 0.0202 0.0140 0.0182 0.0082Table 3
Statistical results of SP for different algorithms.
SP sMOGWO MOGWO MOPSO NSGA-III MOEA/D SPEA2 UF3 Mean 0.5719 1.0807 0.9091 1.0130 1.0233 0.9992 Med 0.53456 1.1129 0.9073 1.0000 1.0196 0.9993 Max 0.7215 1.4434 1.0411 1.0850 1.0747 1.0002 Min 0.4490 0.6002 0.7810 0.9991 1.0000 0.9971  UF4 Mean 0.6098 0.8207 0.7923 0.9991 1.2974 1.5026 Med 0.6217 0.8055 0.7707 0.9940 1.2924 1.4754 Max 0.6551 1.0963 0.9160 1.0616 1.4337 1.9870 Min 0.5428 0.6724 0.6908 0.9464 1.1553 1.1513  UF7 Mean 0.6646 0.8731 0.8649 1.0961 1.0697 0.9995 Med 0.6646 0.8592 0.8386 1.0850 1.0015 0.9996 Max 0.7011 1.2232 1.0611 1.2400 1.5393 1.0760 Min 0.6186 0.5954 0.7292 0.9997 1.0000 0.9553Figure 5
The boxplots of GD value of three algorithms.Figure 6
The boxplots of SP value of three algorithms.The best optimal results of the algorithms in repeated experiments are shown in Figures7, 8, and 9. The comparison is more intuitive.Figure 7
The best optimal results of the different algorithms in repeated experiments on UF3.Figure 8
The best optimal results of the different algorithms in repeated experiments on UF4.Figure 9
The best optimal results of the different algorithms in repeated experiments on UF7.For UF3 test problem, the sMOGWO, NSGA-III, and SPEA2 are the first class in GD value which means the solutions of these methods are closer to the real Pareto front. The GD values of MOPSO, MOGWO, and MOEA/D are high and fluctuating. The sMOGWO method has a better SP value which means the solutions are well distributed. Although the NSGA-III and SPEA2 are closer to the real Pareto front, the SP values are high, which means the solutions of these methods tend to fall into local optimum. We can discover that the NSGA-III and SPEA2 methods fall into local optimum on UF3 test problem intuitively according to Figure7. The solutions of the sMOGWO method are the best.For UF4 test problem, the sMOGWO has the best GD value. The MOPSO, MOGWO, and MOEA/D are in the second class in GD value. The sMOGWO, MOGWO, and MOPSO are in the first class in SP value. From Figure8, the best optimal results of sMOGWO, MOPSO, MOGWO, and MOEA/D are all close to the real Pareto front and are well distributed.For UF7 test problem, the GD values of sMOGWO, MOGWO, NSGA-III, and SPEA2 are in the first class, but the SP value of sMOGWO is better than those of other methods. From Figure9, the best optimal result of all the methods is close to the real Pareto front, but the solutions of MOPSO, MOGWO, MOEA/D, NSGA-III, and SPEA2 are not well distributed. Some methods fall into local optimum. The solutions of the sMOGWO method are the best.Accordingly, the proposed sMOGWO algorithm has good stability under all test functions, and the proposed algorithm performs well in repeated experiments with few poor results.


## 4.1. Test Problems

UF series test problems in CEC 2009 are multimode test functions [29], and there are a lot of locally optimal solutions, so it is suitable for test of multiobjective optimal algorithms. Three typical test problems are chosen as follows: UF3 whose Pareto boundary is concave shape, UF4 whose Pareto boundary is convex shape, UF7 whose Pareto boundary is line. The three UF problems are as follows:UF3 problem,n=30, the search space is 0,1n:(16)minf1=x1+2J14∑j∈J1yj2-2∏j∈J1cos⁡20yiπj+2f2=1-x1+2J14∑j∈J1yj2-2∏j∈J1cos⁡20yiπj+2subjecttoyi=xj-x10.51+3j-2/n-2,j=2,…,n,J1=j∣jisoddand2≤j≤nJ2=j∣jisevenand2≤j≤nUF4 problem,n=30, the search space is [0,1]×-2,2n-1:(17)minf1=x1+2J1∑j∈J1hyif2=1-x12+2J2∑j∈J1hyisubjecttoyi=xj-sin⁡6πx1+jπn,j=2,…,n,ht=t1+e2tJ1=j∣jisoddand2≤j≤nJ2=j∣jisevenand2≤j≤nUF7 problem,n=30, the search space is [0,1]×-1,1n-1:(18)minf1=x15+2J1∑j∈J1yi2f2=1-x15+2J2∑j∈J1yi2subjecttoyi=xj-sin⁡6πx1+jπn,j=2,…,n,J1=j∣jisoddand2≤j≤nJ2=j∣jisevenand2≤j≤n.


## 4.2. Parameters Setting

The principles and the operation modes of the algorithms are different. In order to make the algorithm contrast, the same maximum number of iterations, population size, and Archive size are set up. The parameters of different multiobjective optimization algorithms are shown in Table1.Table 1
The parameters of different multiobjective optimization algorithms.
Description sMOGWO MOGWO MOPSO NSGA-III MOEA/D SPEA2 Maximum number of iterations 1000 1000 1000 1000 1000 1000  Population size nPop=100 nPop=100 nPop=100 nPop=100 nPop=100 nPop=100  Archive size 100 100 100 100 100 100  Inflation rate α=0.1 α=0.1 α=0.1 / / /  Number of grids per dimension nGrid=10 nGrid=10 nGrid=10 / / /  Leader selection pressure β=4 β=4 β=4 / / /  Deletion selection pressure γ=2 γ=2 γ=2 / / /  Inertia weight / / ω=0.5 / / /  Inertia weight damping rate / / ω d a m p=0.99 / / /  Mutation rate / / mu=0.1 mu=0.1 mu=0.1 mu=0.1  Crossover rate / / / 0.5 0.5 0.5  Number of Neighbors / / / / 10 /


## 4.3. Performance Evaluation Indexes

The Generational Distance (GD) and the Spacing (SP) [30] have been used to evaluation the performance of each algorithm. (19)GD=1N∑i=1NDi2whereDi is the Euclidean distance between theith nondominant solution in Pareto solution set and the closest nondominant solution on the real Pareto front.N is the number of Pareto optimal solutions. The smaller value of GD indicates that the Pareto optimal solution is closer to the real Pareto front.(20)SP=df+dl+∑i=1N-1di-d-df+dl+N-1·d-wheredi is the Euclidean distance between the adjacent points of Pareto optimal solution.df anddl are the Euclidean distances between the endpoint of Pareto optimal solution and real Pareto front.N is the number of Pareto optimal solutions. The smaller value of SP indicates that Pareto optimal solution distribution is more homogenized.


## 4.4. Result Analysis

In order to eliminate the contingency, each algorithm is run 30 times independently. The statistical results of the evaluation indexes are shown in Tables2 and 3. In order to reflect the advantages and disadvantages of the optimization results more directly, the boxplots of GD value and SP value of each algorithm are given as in Figures 5 and 6.Table 2
Statistical results of GD for different algorithms.
GD sMOGWO MOGWO MOPSO NSGA-III MOEA/D SPEA2 UF3 Mean 0.0210 0.0542 0.0853 0.0268 0.0561 0.0254 Med 0.0211 0.0540 0.0737 0.0255 0.0571 0.0252 Max 0.0253 0.0900 0.2322 0.0716 0.0829 0.0342 Min 0.0180 0.0315 0.0425 0.0166 0.0299 0.0116  UF4 Mean 0.0239 0.0331 0.0361 0.0516 0.0409 0.0485 Med 0.0241 0.0330 0.0363 0.0521 0.0439 0.0475 Max 0.0257 0.0394 0.0406 0.0622 0.0526 0.0618 Min 0.0234 0.0259 0.0324 0.0415 0.0320 0.0340  UF7 Mean 0.0078 0.0267 0.0468 0.0319 0.0802 0.0313 Med 0.0077 0.0278 0.0424 0.0316 0.0584 0.0301 Max 0.0091 0.1012 0.0881 0.0726 0.1751 0.0608 Min 0.0072 0.0148 0.0202 0.0140 0.0182 0.0082Table 3
Statistical results of SP for different algorithms.
SP sMOGWO MOGWO MOPSO NSGA-III MOEA/D SPEA2 UF3 Mean 0.5719 1.0807 0.9091 1.0130 1.0233 0.9992 Med 0.53456 1.1129 0.9073 1.0000 1.0196 0.9993 Max 0.7215 1.4434 1.0411 1.0850 1.0747 1.0002 Min 0.4490 0.6002 0.7810 0.9991 1.0000 0.9971  UF4 Mean 0.6098 0.8207 0.7923 0.9991 1.2974 1.5026 Med 0.6217 0.8055 0.7707 0.9940 1.2924 1.4754 Max 0.6551 1.0963 0.9160 1.0616 1.4337 1.9870 Min 0.5428 0.6724 0.6908 0.9464 1.1553 1.1513  UF7 Mean 0.6646 0.8731 0.8649 1.0961 1.0697 0.9995 Med 0.6646 0.8592 0.8386 1.0850 1.0015 0.9996 Max 0.7011 1.2232 1.0611 1.2400 1.5393 1.0760 Min 0.6186 0.5954 0.7292 0.9997 1.0000 0.9553Figure 5
The boxplots of GD value of three algorithms.Figure 6
The boxplots of SP value of three algorithms.The best optimal results of the algorithms in repeated experiments are shown in Figures7, 8, and 9. The comparison is more intuitive.Figure 7
The best optimal results of the different algorithms in repeated experiments on UF3.Figure 8
The best optimal results of the different algorithms in repeated experiments on UF4.Figure 9
The best optimal results of the different algorithms in repeated experiments on UF7.For UF3 test problem, the sMOGWO, NSGA-III, and SPEA2 are the first class in GD value which means the solutions of these methods are closer to the real Pareto front. The GD values of MOPSO, MOGWO, and MOEA/D are high and fluctuating. The sMOGWO method has a better SP value which means the solutions are well distributed. Although the NSGA-III and SPEA2 are closer to the real Pareto front, the SP values are high, which means the solutions of these methods tend to fall into local optimum. We can discover that the NSGA-III and SPEA2 methods fall into local optimum on UF3 test problem intuitively according to Figure7. The solutions of the sMOGWO method are the best.For UF4 test problem, the sMOGWO has the best GD value. The MOPSO, MOGWO, and MOEA/D are in the second class in GD value. The sMOGWO, MOGWO, and MOPSO are in the first class in SP value. From Figure8, the best optimal results of sMOGWO, MOPSO, MOGWO, and MOEA/D are all close to the real Pareto front and are well distributed.For UF7 test problem, the GD values of sMOGWO, MOGWO, NSGA-III, and SPEA2 are in the first class, but the SP value of sMOGWO is better than those of other methods. From Figure9, the best optimal result of all the methods is close to the real Pareto front, but the solutions of MOPSO, MOGWO, MOEA/D, NSGA-III, and SPEA2 are not well distributed. Some methods fall into local optimum. The solutions of the sMOGWO method are the best.Accordingly, the proposed sMOGWO algorithm has good stability under all test functions, and the proposed algorithm performs well in repeated experiments with few poor results.


## 5. Case Study

In this section, a case study is presented to show the effect of the proposed algorithm on optimal control of HTGS under multiple operation conditions.The no-load and on-load operation conditions are considered. The model of HTGS and its parameters are provided in Section2. The multiobjective functions refer to (5) and (6). The simulation time is 20 seconds. The sMOGWO, MOGWO, MOPSO, MOEA/D, NSGA-III, and SPEA2 algorithms are utilized for optimal control. The maximum number of iterations is 100. The other parameters of the three algorithms are set as in Table 1. Each algorithm is run 30 times independently. The best optimal results of each algorithms in repeated experiments are shown in Figure 10.Figure 10
The Pareto front of different algorithms.The solutions of MOPSO method are far away from Pareto front. The solutions of sMOGWO, MOGWO, MOPSO, MOEA/D, NSGA-III, and SPEA2 are similar to the Pareto front, but the solutions of sMOGWO are more evenly distributed than those of other methods.In order to analyze the detailed control transient process of the optimization results, some representative control strategies are chosen from the Pareto optimal solution set as shown in Table4. It can be seen that the optimization of ITAE under no-load and ITAE under on-load are opposites. Transient process of some of the control strategies is shown in Figure 11 to exhibit control effects more intuitively.Table 4
Some representative control strategies are chosen from the Pareto optimal solution set.
Control strategy Parameters of HGTS Performance indexes K p K i K d ITAE under no-load ITAE under on-load 1 8.750 5.173 1.921 0.2352 0.0176  2 9.352 4.612 1.602 0.2156 0.0220  3 9.791 3.864 1.216 0.1925 0.0311  4 10.171 3.514 1.028 0.1752 0.0402  5 10.473 3.191 0.759 0.1585 0.0517  6 10.967 2.595 0.271 0.1284 0.0823Figure 11
Transient process by different control strategies.All the three control strategies have good control stability under multiple operation conditions. Strategy 1 is most stable under on-load operation condition, but it has bigger overshoot than the other two strategies under no-load operation condition. Strategy 3 has the smallest overshoot among the three strategies under no-load operation condition, but its stability is the worst among the control strategies under on-load operation condition. Strategy 2 is a compromise control strategy.Thus, a Pareto solution set can be got after one optimization which can be suitable for multiple operation conditions. And the solutions have different emphasis on objective functions. According to the specific requirements of the HTGS, the decision maker can select one or some satisfactory solutions from this solution set as the final solution, so that the HTGS can obtain better control quality.


## 6. Conclusion

The control strategies of hydraulic turbine governing system need to consider the multiple operation conditions. A multiobjective optimal function under different operation conditions is proposed in this paper to solve this control problems.In order to optimize the multiobjective problems more effectively, a novel MOGWO algorithm based on searching factor (sMOGWO) is proposed. The sMOGWO method is verified with several UF test problems compared to MOGWO, MOPSO, MOEA/D, NSGA-III, and SPEA2. The sMOGWO provided better solution compared with its competitors. And the proposed algorithm has good stability under all test functions, and the proposed algorithm performs well in repeated experiments with few poor results.A case study has been designed to test the control quality of the control strategies which are got by the proposed method. The experimental results have confirmed that the control strategies perform well under multiple operation conditions.



---

*Source: 3745924-2019-05-28.xml*