# A Hybrid Multiobjective Optimization Based on Nondominated Sorting and Crowding Distance, with Applications to Wave Energy Converters

Authors: John Saveca; Yanxia Sun; Zenghui Wang
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

Multiobjective evolutionary algorithm based on decomposition (MOEA/D) has attracted a lot of attention since it can handle multiobjective problems (MOP) with a complicated Pareto front. The procedure involves decomposing a MOP into single subproblems, which are eventually optimized simultaneously based on the MOP neighborhood information. However, the MOEA/D strategy tends to produce a distributed optimization that is not of good quality in some problems with complex Pareto optimal front, such as problems with a long tail and sharp peak, common in real-world situations. This paper proposes an improved MOEA/D to enhance the distributed optimization quality and minimize its complexity while accelerating the optimization to get a better solution. The improved method is achieved by incorporating a Hybrid Differential Evolution/Particle Swarm Optimization algorithm and a hybrid operator based on nondominated sorting and crowding distance algorithm. This incorporation takes place in the mutation generator and initial population part of the original MOEA/D algorithm. Simulations and comparisons are carried out based on some MOP benchmark functions to verify the proposed method’s performance. The experimental results show that the proposed method achieves better performance compared to other algorithms. Furthermore, the proposed method is also applied to optimize the multiobjective wave energy converter model to maximize power per year and minimize cost per unit power.


---

## Body

## 1. Introduction

The study of MOEA/D has been a hot topic lately [1]. It has proved to be amongst the highly competitive and dominating Multiobjective Evolutionary Algorithms (MOEAs) alongside other algorithms such as Multiobjective Differential Evolution (MODE) [2], Multiobjective Particle Swarm Optimization (MOPSO) [3], and Nondominated Sorting Genetic Algorithm-II (NSGA-II) [4]. Like other evolutionary algorithms (EAs), MOEA/D is also a population-based MOEA, and it is inspired by a decomposition mechanism that decomposes a multiobjective problem (MOP) into a set of single-objective optimization subproblems that are equivalent to the population of candidate solutions (on other MOEAs), where each subproblem is assigned with a unique weight vector [1]. The subproblems are then optimized simultaneously in a cooperative manner [5]. Each subproblem’s optimal solution becomes a Pareto optimal solution to the MOP being optimized [6].There have been some evolutionary algorithms (EAs) developed recently by several researchers. Monarch Butterfly Optimization (MBO) is one of them, and it is a class of swarm intelligence algorithms inspired by the migration behavior of monarch butterflies [7]. Another recent EA is Slime Mould Algorithm (SMA), proposed by Li et al. [8] to solve single-objective optimization problems, known for its efficient global search ability [9]. Harris hawks optimization (HHO) is a nature-inspired metaheuristic recently proposed by [10]. The main inspiration of the algorithm is the cooperative behavior, and chasing style of Harris’ hawks in nature is called surprise pounce, according to [10]. However, most of these recently developed metaheuristics are still being employed for single-objective optimization. However, this study focuses on multiobjective optimization problems; at the moment, only SMA has been developed to handle multiobjective problems in the form of Multiobjective Slime Mould Algorithm [11] and it is also used to compare its performance with the proposed algorithm.There are several studies developed based on the hybrid DE/PSO algorithm. Many of these studies focused on optimizing a single objective instead of multiple objectives, covered in [12–18]. However, this study aims to develop an improved version of MOEA/D based on hybrid operators of DE, PSO, and CD with NdS.Many researchers have praised the MOEA/D since its introduction by Zhang and Li in 2007 due to its advantages [19]. According to [19], MOEA/D provides a path that presents the decomposition techniques into multiobjective evolutionary optimization that is easy to understand and efficient. Since MOEA/D optimizes many single-objective optimization problems simultaneously in parallel rather than the whole MOP, issues such as the assignment of fitness and maintenance of diversity that cause difficulties for other MOEAs could be handled easily in the MOEA/D framework [19]. Besides having many advantages compared to other MOEAs, MOEA/D has its disadvantages. It has been highlighted by [5, 20] that the MOEA/D does not produce good quality and a well-distributed optimization in problems with complicated Pareto optimal front shapes. Such problems are problems with long tail or sharp peak [5, 20], which qualifies the complexity of real-world problems, like the Sharp Edge FIR Filters Design problem [21], Mitigation of Power Oscillations problem [22], and a Long Tail problem in Recommender System [23]. Several researchers proposed improvements on the MOEA/D to overcome the issues mentioned above. He et al. proposed the normalization of objective space technique in an integrated manner where two normalization methods (MOEA/D-2N) are employed simultaneously to remedy the performance diversity and deterioration of MOEA/D when applied to the multiobjective problems with real-world characteristics [24]. Li and Chen [25] proposed the MOEA/D improvement by adaptive weight vector and matching strategy (MOEA/D-AVM) to overcome the decline of the population diversity and the performance degradation of the solution set when solving a MOP of discontinuous Pareto front [25]. Wang et al. proposed using two reference points, namely, the ideal point and the nadir point, to impact the MOEA/D algorithm performance during optimization of MOPs’ with complex Pareto optimal front [26]. For this study, incorporating a Hybrid Differential Evolution/Particle Swarm Optimization algorithm and a hybrid operator based on nondominated sorting and crowding distance algorithm into MOEA/D algorithm can improve the optimal solution, optimization quality, and distribution on MOPs’ with complex Pareto front.Renewable energies are getting popular nowadays as the source that will sustain the future generation of electric power. With continuous air pollution and forever increments of per-unit energy costs of fossil fuel-generated energy, the wave energy converter (WEC) devices are determined to present clean, safe, and cheap energy to society [27]. The actual application of the WEC device is to convert the wave energy through its motion into electricity, based on various working principles that classify WECs into different categories [27]. The categories of WEC devices are point absorbers, attenuators, and terminators. According to [28] by the year 2009, more than 1000 WEC devices were already researched and developed. These stats indicate how much of a research magnet the WEC topic is to the industries and research institutions. This gives an indication that the study of WECS still has some difficulties that are yet to be tackled, such as energy conversion efficiency, high maintenance costs, and installation costs as argued by [29]. High maintenance costs are due to extreme conditions that influence harsh waves, leading to the damaging of WEC devices [30]. Another shortcoming of WEC devices is the irregularity of ocean waves which produces fluctuating output power. Therefore, some means of power compensating systems are needed in place to achieve a smooth power output [30]. These difficulties come as a result of the WEC study involving multidisciplinary components for harvesting power from ocean waves [28]. For instance, the WEC performance is highly determined by the PTO (power take-off) system, the hydrodynamic design, and dynamics and control in an attempt to increase the WEC efficiency for different sea states. This multidisciplinary nature results in optimization difficulties of WEC devices [28]. According to [31], the annual total average of wave energy calculated at a water depth of 60 m of the U.S. coastlines is estimated at 2100 Terawatt hours (TWh). WEC Energy has the highest average power density of 25 KW/m2 compared to other renewable sources such as solar and wind energy, which both have a low average power density of 1 kW/m2 each, at peak solar insolation and 12 m/s, respectively [32]. This makes the WEC system one of the reliable power harvesters for a more extended period should the current fossil fuels cease to exist.As mentioned above, WECs come in different types, the point absorber type offers many benefits compared to other types, and it has been researched more than other types. The point absorber benefits are the low complexity and high reliability, which are significant for an offshore platform, as stated by [30]. This paper proposes to optimize the point absorber WEC device with the objectives of maximizing the annual generated power and minimizing the cost of the per-unit energy generated by applying an improved version of MOEA/D. During the optimization of the abovementioned WEC problem, the point absorber heaving buoy shape, size, and submergence will be taken into account. The constraints of designing the point absorber WEC device for this research are the spherical shape buoy radius variable boundaries and the cylindrical buoy shape length and radius variables boundaries [32]. Various researchers have conducted many studies on point absorber optimization. A parameter study and optimization of two-body wave energy converters was conducted by [28]. Their research presented a parametric study, design, and optimization procedures that studied the effects of the different parameters [28]. Another study on point absorbers was conducted by [31]. Their study investigated a novel design of the WEC harvester that utilized a flywheel energy storage system to produce an increased power generation. An optimization of heaving buoy WEC based on a combined numerical model was proposed by [27]. The study focused on the mass, dimension, hydrodynamic response, and power take-off (PTO) damping for heaving buoy wave energy converters (WEC). For this study, the energy extraction will be maximized while keeping the cost of per-unit power at a minimum by applying the proposed HMOONDS-CD algorithm when a point absorber WEC is employed.While DE is inspired by the existence of a special type of difference vector [33], and PSO is inspired by the social behavior of birds flocking or fish schooling [34, 35], both algorithms have shown a great track of effective and efficient performance in their history of problem optimization. In this research, an improved MOEA/D algorithm based on Particle Swarm Optimization/Differential Evolution Hybrid with fixtures of nondominated sorting and crowding distance is proposed to overcome the abovementioned challenges. The MOPs with complicated Pareto optimal front will be used to evaluate the performance of the proposed improved algorithm. The improved algorithm will further be applied to the point absorber WEC optimization problem. The rest of the paper is organized as follows. Section 2 discusses the formation of MOP and the background of the traditional MOEA/D and its framework. Section 3 discusses the improved MOEA/D based on hybrid PSO/DE with fixtures of nondominated sorting and the crowding distance.The section also discusses the framework of the improved MOEA/D with its fixture. Section4 discusses the experiment setup and the performance matrix indicators to be used in this research. The section also outlines the MOPs to be used in the experiment with their characteristics. Section 5 discusses the results of simulations of the experiment. Section 6 concludes the research findings. Section 7 gives the future works.


## 2. Multiobjective Optimization Problem and Original MOEA/D Background

### 2.1. Multiobjective Optimization Problem

A minimized multiobjective optimization problem (MOP) is defined as follows:(1)minimizeFx=f1x,f2x,…,fmxTs.t.x∈Ω,f∈Rm,where x=x1,x2,…,xn∈Ω is n-dimensional decision variables vector in the feasible search space, Ω is the feasible search region, and n is the dimensions of the decision vector. Rm is the objective region, and m is the number of objectives [36]. In multiobjective optimization, the objectives in (1) usually conflict with each other when being optimized; therefore, the improvement of one objective results in the deterioration of at least one other objective [19, 37]. This means that there is no point in the decision space (Ω) that can minimize or maximize all the objectives at the same instance. In that case, the MOEA must find Pareto optimal solutions, which are the best trade-off solutions to help decision-makers comprehend the trade-off connection between the contradicting objectives and then make the final decision [1, 19, 37]. Following are Pareto optimal concepts related conditions:Condition 1: Letq,p∈Ω, be the two solutions, q is said to dominate p if, and only if, qi≤pi for every i∈1,…,m and qj<pj for at least one index j∈1,…,m, denoted as q<p [1, 38].Condition 2: A feasible solutionx∗∈Ω of problem (1) is said to be Pareto optimal if there is no point x∈Ω such that Fx dominates Fx∗.Fx∗ and it is called Pareto optimal (objective) vector [19].Condition 3: The set of all the Pareto optimal points is called the Pareto optimal setP∗, denoted as P∗=x∈Ω|Fx,Fx∗.Fx∗ [37].Condition 4: TheP∗ image in the objective region is called the Pareto front PF∗, denoted as PF∗=fx|,x∈P∗ [37].


### 2.2. Multiobjective Evolutionary Algorithm Based on Decomposition Framework

The decomposition mechanism is the driving component in MOEA/D [39]. It dissolves a multiobjective problem into a number of single-objective problems which are optimized simultaneously. Each subproblem is optimized by using information from its several neighboring subproblems to effectively reduce the algorithm’s computational complexity. In this paper, the Tchebycheff decomposition method will be used [20]. The Tchebycheff approach is a modifying factor when the decision maker’s preferred information proceeds as the main algorithm progresses to find a set of nondominated solutions [40]. The preferred information is represented based on user-specified inputs [40]. One of the preferred information is the ideal reference solution (Zref) to the MOP [40]. The other preferred information is the weight vector λi that assigns the relative preference of the user [40]. The objective function of the subproblem to be minimized is defined as follows [19]:(2)minimizegtex|λ,Z∗=maxλjifjx−Zj∗,Subject tox∈Ω,where gte is the subproblem objective function to be minimized, λi is a weight vector, the value of k ranges from 1 to N, where N is the number of populations, and j is the number of single-objective optimization subproblems. The weight vector λi is distributed equally, and the neighborhood size T weight vectors closest to λi in the weight vector are generally selected as neighbors, and the neighbors of the ith subproblem are denoted as Bi [19]. fjx is the jth objective [20].Z∗=Z1∗,…,Zn∗T is the ideal reference solution, i.e., Zj∗=minfjx,x∈Ω for each j=1,…,n.The following is the general framework of MOEA/D according to Zhang and Li [19]. Suppose the Tchebycheff method is selected for the decomposition approach and λ1,…,λN is a group of even distributed weight vectors and Z∗ is the reference solution. The problem on (1) becomes the PF approximation problem that is decomposed into N of single-objective optimization subproblems. The j th subproblem’s objective function is given as follows:(3)gtex|λj,Z∗=max1≤i≤mλjifjx−Zj∗,where λj=λ1j,…,λmjT. In this case, all the N objective functions are optimized simultaneously in a single run, and λ is influential to gte. At this instant, with the scalar objective function gtex|λi,Z∗ the optimal solution is relatively close to the optimal solution of the scalar objective function gtex|λj,Z∗ when their weight vectors λi and λj are close to each other. Therefore, other scalar objective functions (gtex|λN,Z∗ with weight vectors close to that of λi will have immediate optimal solutions that serves as a massive boost to the MOEA/D algorithm. The i th subproblem’s neighborhood is made up of all the subproblems with the weight vectors from the neighborhood ofλi. The population consists of the optimal solution obtained so far for each subproblem. Only the current solutions to its neighboring subproblems are exploited for optimizing a subproblem in MOEA/D [19]. Figure 1 below shows the MOEA/D flowchart according to [5].Figure 1
Flowchart of pseudo for the original MOEA/D algorithm [5].


## 2.1. Multiobjective Optimization Problem

A minimized multiobjective optimization problem (MOP) is defined as follows:(1)minimizeFx=f1x,f2x,…,fmxTs.t.x∈Ω,f∈Rm,where x=x1,x2,…,xn∈Ω is n-dimensional decision variables vector in the feasible search space, Ω is the feasible search region, and n is the dimensions of the decision vector. Rm is the objective region, and m is the number of objectives [36]. In multiobjective optimization, the objectives in (1) usually conflict with each other when being optimized; therefore, the improvement of one objective results in the deterioration of at least one other objective [19, 37]. This means that there is no point in the decision space (Ω) that can minimize or maximize all the objectives at the same instance. In that case, the MOEA must find Pareto optimal solutions, which are the best trade-off solutions to help decision-makers comprehend the trade-off connection between the contradicting objectives and then make the final decision [1, 19, 37]. Following are Pareto optimal concepts related conditions:Condition 1: Letq,p∈Ω, be the two solutions, q is said to dominate p if, and only if, qi≤pi for every i∈1,…,m and qj<pj for at least one index j∈1,…,m, denoted as q<p [1, 38].Condition 2: A feasible solutionx∗∈Ω of problem (1) is said to be Pareto optimal if there is no point x∈Ω such that Fx dominates Fx∗.Fx∗ and it is called Pareto optimal (objective) vector [19].Condition 3: The set of all the Pareto optimal points is called the Pareto optimal setP∗, denoted as P∗=x∈Ω|Fx,Fx∗.Fx∗ [37].Condition 4: TheP∗ image in the objective region is called the Pareto front PF∗, denoted as PF∗=fx|,x∈P∗ [37].


## 2.2. Multiobjective Evolutionary Algorithm Based on Decomposition Framework

The decomposition mechanism is the driving component in MOEA/D [39]. It dissolves a multiobjective problem into a number of single-objective problems which are optimized simultaneously. Each subproblem is optimized by using information from its several neighboring subproblems to effectively reduce the algorithm’s computational complexity. In this paper, the Tchebycheff decomposition method will be used [20]. The Tchebycheff approach is a modifying factor when the decision maker’s preferred information proceeds as the main algorithm progresses to find a set of nondominated solutions [40]. The preferred information is represented based on user-specified inputs [40]. One of the preferred information is the ideal reference solution (Zref) to the MOP [40]. The other preferred information is the weight vector λi that assigns the relative preference of the user [40]. The objective function of the subproblem to be minimized is defined as follows [19]:(2)minimizegtex|λ,Z∗=maxλjifjx−Zj∗,Subject tox∈Ω,where gte is the subproblem objective function to be minimized, λi is a weight vector, the value of k ranges from 1 to N, where N is the number of populations, and j is the number of single-objective optimization subproblems. The weight vector λi is distributed equally, and the neighborhood size T weight vectors closest to λi in the weight vector are generally selected as neighbors, and the neighbors of the ith subproblem are denoted as Bi [19]. fjx is the jth objective [20].Z∗=Z1∗,…,Zn∗T is the ideal reference solution, i.e., Zj∗=minfjx,x∈Ω for each j=1,…,n.The following is the general framework of MOEA/D according to Zhang and Li [19]. Suppose the Tchebycheff method is selected for the decomposition approach and λ1,…,λN is a group of even distributed weight vectors and Z∗ is the reference solution. The problem on (1) becomes the PF approximation problem that is decomposed into N of single-objective optimization subproblems. The j th subproblem’s objective function is given as follows:(3)gtex|λj,Z∗=max1≤i≤mλjifjx−Zj∗,where λj=λ1j,…,λmjT. In this case, all the N objective functions are optimized simultaneously in a single run, and λ is influential to gte. At this instant, with the scalar objective function gtex|λi,Z∗ the optimal solution is relatively close to the optimal solution of the scalar objective function gtex|λj,Z∗ when their weight vectors λi and λj are close to each other. Therefore, other scalar objective functions (gtex|λN,Z∗ with weight vectors close to that of λi will have immediate optimal solutions that serves as a massive boost to the MOEA/D algorithm. The i th subproblem’s neighborhood is made up of all the subproblems with the weight vectors from the neighborhood ofλi. The population consists of the optimal solution obtained so far for each subproblem. Only the current solutions to its neighboring subproblems are exploited for optimizing a subproblem in MOEA/D [19]. Figure 1 below shows the MOEA/D flowchart according to [5].Figure 1
Flowchart of pseudo for the original MOEA/D algorithm [5].


## 3. Proposed MOEA/D Based on Hybrid DE/PSO and Nondominated Sorting with Crowding Distance

Diversity, efficiency, complexity, and optimality are the elements of the algorithm to consider while searching for the optimum solution in MOP and the proposed improved MOEA/D based on hybrid DE/PSO with NdS and CD fixtures will consider those elements. DE mutation phase generates the best new solutions from the population. In contrast, the velocity vector update phase selects the best possible solutions, and it is much effective on optimization acceleration enhancement. Therefore with these strategies made hybrid, the proposed algorithm has the ability to solve the problems mentioned above, with the help of nondominated sorting [41], which uses Pareto dominance to sort a group of solutions for the MOP and the crowding distance, which acts as a metric of the second order from the primary nondominated sorting algorithm. The proposed algorithm will handle the constraints on the benchmark functions by limiting the function variable value to the nearest border when the variable is out of range. For the improved MOEA/D to have a better selection of the nondominated solutions from the initial population solutions, the NdS fixture together with the CD is used to select the best initial solutions from the population and to select the best fitness initial value for updating the reference point Z∗, the nadir point [5] Znad, and the initialization of the subproblems.


### 3.1. The Nondominated Sorting and Crowding Distance Application on the Proposed HMOONDS-CD

In the improved MOEA/D, for each solution on the populationP space Ω and the Fitness Value FV space Ω , the NdS will compute the number of all the solutions Np which dominate the solution P and the number of all the solutions NFV that dominate the solution FV. NdS will also compute the groups of solutions Gp and GFV that the solution P and solution FV dominate, respectively [2, 42]. The computation takes place simultaneously and in parallel for solution P and solution FV. All the first rank-dominated front solutions for P and FV will be set to a count zero (0). Therefore, for each solution P with Np=0 and each solution FV with NFV=0, their members qp and qFV of their groups Gp and GFV are visited, respectively, and their dominations are reduced by count one (1). At this stage, if the domination count for any member becomes zero, then that member is moved to a disjoint list denoted by Q. At this moment, list Q is the second rank-dominated front. The same procedure is continued for identifying the third member and until all the members are identified. The CD will approximate the density of the solutions from the NdS in P∈Ω and FV∈Ω. This will be achieved by computing the arithmetic mean distance value between two solutions on either side of the individual target solution along each of the objectives [41, 43].


### 3.2. Differential Evolution/Particle Swarm Optimization Hybrid

DE and PSO are the most preferred evolutionary algorithms [44, 45], and they have been transformed into a robust state-of-the-art algorithm to date. This is due to their efficiency and effectiveness in application and implementation compared to other mathematical and evolutionary algorithms [34, 46]. With the abovementioned advantages both algorithms possess, they are made a hybrid operation for the improvement of the HMOONDS-CD based on the abovementioned nondominated sorting and crowding distance algorithms to be applied on WEC optimization. The DE-PSO hybrid algorithm has been used by several researchers, such as [47], for the suspended sediment load estimation application, where the DE-PSO hybrid was coupled with multilayer perceptron algorithm. The DE-PSO hybrid has also been applied to design and implement sharp edge FIR filters [21] and to mitigate power oscillations based on SSSC damping controller [22]. The mutation and crossover phases are the driving forces in the success of DE efficient operation, while the velocity vector update phase and the position vector phase are the engines of the PSO convergence speed and effectiveness during operation. During the DE phase, three individuals r1,r2,r3r1≠r2≠r3≠i are randomly selected from a set of subproblems Bi in order to generate the new solution y. DE/Rand/1 mutation strategy is selected, and it is expressed as follows:(4)yi,=Xr1,+FiXr2,−Xr3,,where r1,r2,r3 are three individuals randomly chosen from the set Bi while Xr2,—Xr3, are the different vectors to mutate the parent and Xr1 is the current mutant vector [46]. During the mutation search session, the crossover phase complements the DE/Rand/1 mutation strategy by exchange information between different individuals in the current search session, expressed as follows:(5)yj=Xj,i,g......frandj0,1≤CRiorj=jrand,Xj,i,g.........otherwise,where CR ∈ [0, 1] is the crossover probability, a user-defined value that controls the parameter values fraction that are copied from the mutant [48]. The new solution obtained from DE mutation is further updated based on the Swarm Velocity and position update on the PSO phase. This strategy enhances the algorithm optimization speed for fast selection of the best possible solutions to be fed to the decomposition stage. The PSO particle best position pBest is set from the DE new solution y while the PSO index is set as a maximum function of the optimization problem with respect to the DE new solution y in order to establish the PSO group best gBest parameter. The new updated y¯ from PSO operation is established based on the velocity swarm update where c1 and c2 are learning factor 1 and global learning factor 2, respectively [47]. Parameters r1 and r2 are random parameters with r1 ranging between [pBest,y] values while r2 ranges between [gBest,y] values. Parameter ω is the inertia weight, and v is the population-based velocity by PSO independent population size and MOP problem boundaries. The PSO new solution y¯ is further updated based on PSO position vector as expressed in (10). Figure 2 below shows the proposed HMOONDS-CD flowchart.(6)pBest=y,(7)index=maxfjx,y,(8)gBest=yindex,(9)y¯it+1=ωvit+c1r1pBestit−yit+c2r2gBestit−yit,(10)yiPSO=y¯it+vit.Figure 2
Flowchart for pseudo of the improved HMOONDS-CD algorithm.


## 3.1. The Nondominated Sorting and Crowding Distance Application on the Proposed HMOONDS-CD

In the improved MOEA/D, for each solution on the populationP space Ω and the Fitness Value FV space Ω , the NdS will compute the number of all the solutions Np which dominate the solution P and the number of all the solutions NFV that dominate the solution FV. NdS will also compute the groups of solutions Gp and GFV that the solution P and solution FV dominate, respectively [2, 42]. The computation takes place simultaneously and in parallel for solution P and solution FV. All the first rank-dominated front solutions for P and FV will be set to a count zero (0). Therefore, for each solution P with Np=0 and each solution FV with NFV=0, their members qp and qFV of their groups Gp and GFV are visited, respectively, and their dominations are reduced by count one (1). At this stage, if the domination count for any member becomes zero, then that member is moved to a disjoint list denoted by Q. At this moment, list Q is the second rank-dominated front. The same procedure is continued for identifying the third member and until all the members are identified. The CD will approximate the density of the solutions from the NdS in P∈Ω and FV∈Ω. This will be achieved by computing the arithmetic mean distance value between two solutions on either side of the individual target solution along each of the objectives [41, 43].


## 3.2. Differential Evolution/Particle Swarm Optimization Hybrid

DE and PSO are the most preferred evolutionary algorithms [44, 45], and they have been transformed into a robust state-of-the-art algorithm to date. This is due to their efficiency and effectiveness in application and implementation compared to other mathematical and evolutionary algorithms [34, 46]. With the abovementioned advantages both algorithms possess, they are made a hybrid operation for the improvement of the HMOONDS-CD based on the abovementioned nondominated sorting and crowding distance algorithms to be applied on WEC optimization. The DE-PSO hybrid algorithm has been used by several researchers, such as [47], for the suspended sediment load estimation application, where the DE-PSO hybrid was coupled with multilayer perceptron algorithm. The DE-PSO hybrid has also been applied to design and implement sharp edge FIR filters [21] and to mitigate power oscillations based on SSSC damping controller [22]. The mutation and crossover phases are the driving forces in the success of DE efficient operation, while the velocity vector update phase and the position vector phase are the engines of the PSO convergence speed and effectiveness during operation. During the DE phase, three individuals r1,r2,r3r1≠r2≠r3≠i are randomly selected from a set of subproblems Bi in order to generate the new solution y. DE/Rand/1 mutation strategy is selected, and it is expressed as follows:(4)yi,=Xr1,+FiXr2,−Xr3,,where r1,r2,r3 are three individuals randomly chosen from the set Bi while Xr2,—Xr3, are the different vectors to mutate the parent and Xr1 is the current mutant vector [46]. During the mutation search session, the crossover phase complements the DE/Rand/1 mutation strategy by exchange information between different individuals in the current search session, expressed as follows:(5)yj=Xj,i,g......frandj0,1≤CRiorj=jrand,Xj,i,g.........otherwise,where CR ∈ [0, 1] is the crossover probability, a user-defined value that controls the parameter values fraction that are copied from the mutant [48]. The new solution obtained from DE mutation is further updated based on the Swarm Velocity and position update on the PSO phase. This strategy enhances the algorithm optimization speed for fast selection of the best possible solutions to be fed to the decomposition stage. The PSO particle best position pBest is set from the DE new solution y while the PSO index is set as a maximum function of the optimization problem with respect to the DE new solution y in order to establish the PSO group best gBest parameter. The new updated y¯ from PSO operation is established based on the velocity swarm update where c1 and c2 are learning factor 1 and global learning factor 2, respectively [47]. Parameters r1 and r2 are random parameters with r1 ranging between [pBest,y] values while r2 ranges between [gBest,y] values. Parameter ω is the inertia weight, and v is the population-based velocity by PSO independent population size and MOP problem boundaries. The PSO new solution y¯ is further updated based on PSO position vector as expressed in (10). Figure 2 below shows the proposed HMOONDS-CD flowchart.(6)pBest=y,(7)index=maxfjx,y,(8)gBest=yindex,(9)y¯it+1=ωvit+c1r1pBestit−yit+c2r2gBestit−yit,(10)yiPSO=y¯it+vit.Figure 2
Flowchart for pseudo of the improved HMOONDS-CD algorithm.


## 4. Heaving Buoy Pont Absorber WEC Optimization Problem

Ocean wave energy is one of the world’s most robust forms of energy, and extracting it from oceans has proven to be very strenuous [30]. The high cost of wave energy limits it to be competitive on the market, and this limits big companies to participate in this sector [30, 49]. However, if the research and development of this form of energy is thoroughly done to the point where the high cost of the energy is minimized and the power is efficiently extracted, more interest from big industries would be achieved. This section focuses on optimizing the point absorber to maximize the annual power generated and minimize the Levelized Cost of Energy (LCOE), taking buoy shape and radius height into account. HMOONDS-CD will again be applied to this problem. Following is the expression of the mathematical model of the wave energy converter system to be optimized.The buoy motion is presented by Newton’s second law of motion as expressed on (11), where m is the comprehensive buoy mass, and a is the acceleration of the buoy, and F⇀buoy is the total sum of the forces that act on the buoy [50].(11)mbuoy⋅a⇀=∑F⇀buoy.In this case, the forces acting on the buoy are identified as the buoy weight, the buoyancy force, the fluid drag force, the wave excitation force, and the cable tension force, and they are illustrated in Figure3 [49, 50].Figure 3
Forces acting on the heaving buoy [49].


### 4.1. Buoy Weight

The weight being the simplest force that acts on the buoy, it is expressed on (12), where mbuoy is the mass of the entire buoy system, and g is the gravitational constant [49, 51].(12)W=mbuoy⋅g.The buoy mass of the point absorber WEC system can be increased by attaching the weights if needed in the instances where the buoy is not heavy enough and skips out of the water when it impacts with large waves [49, 51]. It should also be noted that the weight should not exceed the buoyancy force of a fully submerged buoy, or else it would cause the buoy to drop down the ocean [49, 51].


### 4.2. Buoyancy Force

The buoyancy force is exerted by the fluid in the upward direction of a partially or fully submerged object, opposing the object’s weight. The phenomenon is influenced by the pressure difference between the bottom and the top of the submerged object that acts on its surfaces in a fluid [52, 53], and it is expressed on (13), where Vsubmerged is the buoy submerged portion volume in the sea and ρw is the water density [49, 54].(13)Fbuoyancy=ρw⋅g⋅Vsubmerged.In this case, the submerged volume change is influenced by the buoy position change. The volume of the submerged portion of the buoy is expressed as follows:(14)Vsubmerged=π⋅R2⋅Zwave−z,zwave−h≤z≤zwave,where h and R are cylindrical height and radius, respectively, Zwave is the wave height, and z coordinate represents the buoy depth from the perspective of the bottom surface [53, 54]. If z=0, the buoy bottom surface is level to ocean top surface; if z=−6, the buoy is submerged 6 meters below the ocean surface. In the case of spherical shape, the submerged volume is expressed in(15)Vsubmerged=43⋅π⋅R3⋅Zwave−z.The volume of the cylindrical buoy is determined by its radius and height, while the spherical shaped buoy volume is only determined by its radius. These geometric construction characteristics can have a different influence on the energy extraction of the WEC system. With that factor, the experiment will determine which buoy shape has the upper hand on wave energy extraction.


### 4.3. Fluid Drag Force

Drag force becomes present when the body and a fluid interact to produce the relative motion, and it acts opposite to body motion direction [49, 52]. The selected expression of the drag force is based on the Reynolds number [44], Re, which computes the flow condition, as expressed on (16):(16)Re=ρw⋅v⋅Lμw,where μw is the dynamic water viscosity, v is the relative buoy velocity from the water perspective, and L is the characteristic linear dimension [55]. Drag force is expressed as follows:(17)Fdrag=12ρw⋅v2⋅Cd⋅A,where Cd is the drag coefficient and A is the surface area influenced by the drag force [48, 50].


### 4.4. Excitation Force

The most intricate force exerted on the buoy is the wave excitation force. It is a superposition of the wave’s incident, diffraction, and radiation potentials [56]. If the excitation force is assumed to be equivalent to the incident potential, it is referred to as the Froude-Krylov force [39, 56], expressed on (18).(18)FFK=−ρw∬sub∂ΦI∂tn^⋅dS.The partial derivative from the perspective of time is done for the potential of the incident waveΦI, while the surface integral is taken over the submerged surfaces, where n^ is a normal vector unit to the surface, and dS is the dot product with the incremental surface area [49, 52]. The incident potential is expressed on (19), where A is the incident wave amplitude, ωw is the angular frequency in radian per second, k is the angular wavenumber in radian per meter, z is the depth of the submerged buoy, t is time, x is the horizontal position, and Re is the operator that keeps the actual component of a complex number [39].(19)ΦI=A⋅ωwkekzReieiωwt−kx.After the Froude-Krylov formula is simplified, the derived excitation force is estimated as in (20), where vwave is the vertical velocity of the wave [39].(20)FEXC=ρw⋅πR2k⋅vwave⋅k⋅z.


### 4.5. Mechanical System Model

The model of a mechanical system combines all the forces that act on the cable-driven pulley. Newton’s second law of motion is used to express the rotational motion, as expressed in (21), where I is the moment of inertia of the whole system being rotated by the pulley, α is the pulley angular acceleration, and τpulley are the torques that are laid on the pulley [49, 50].(21)I⋅a=∑τpulley.Equation (22) shows the resistive torque τres, where βback and βfric are the coefficients used for back-torque and mechanical friction, respectively, and CLoad is the coefficient of the load [49].(22)τres=CLoad⋅βback+βfric⋅ω.


### 4.6. Power Output and Per-Unit Cost

The WEC system output power from the linear generator is approximated to the square of the rotor’s angular velocity [49]. The instantaneous power P, generated by the generator, is given by (23), where Cpower is the generator-dependent power coefficient [49].(23)P=CLoad⋅Cpower⋅ω2.Objective one is derived based on (23), where fx1 is equated to the sum of the power generated by the liar generator. The WEC system’s annual generated energy is expressed on (24). During the multiobjective optimization of the WEC system, objective one will be maximized while objective 2 as expressed on (27) will be minimized.(24)Objective1:fx1=∑Pave,were Pave is the average power generated by a linear generator that is expected to be maximized.


### 4.7. Annual Energy

(25)Eannual=P⋅Hoursdaily⋅daysyearly.LCOE, expressed on (26), is a measure of the average net present cost of energy for a generating plant lifetime [54]. In this case, the generating period is set to 1 year to determine the average cost of the per-unit energy vs. the total annual power generated by the optimized WEC system.(26)LCOE=Eannualm−0.5,(27)Objective2:fx2=LCOE,where LCOE is the Levelized Cost of Energy to be minimized, Eannual is the total annual energy, m is the buoy mass in kilograms, and −0.5 is inserted for the optimization process based on [54].


## 4.1. Buoy Weight

The weight being the simplest force that acts on the buoy, it is expressed on (12), where mbuoy is the mass of the entire buoy system, and g is the gravitational constant [49, 51].(12)W=mbuoy⋅g.The buoy mass of the point absorber WEC system can be increased by attaching the weights if needed in the instances where the buoy is not heavy enough and skips out of the water when it impacts with large waves [49, 51]. It should also be noted that the weight should not exceed the buoyancy force of a fully submerged buoy, or else it would cause the buoy to drop down the ocean [49, 51].


## 4.2. Buoyancy Force

The buoyancy force is exerted by the fluid in the upward direction of a partially or fully submerged object, opposing the object’s weight. The phenomenon is influenced by the pressure difference between the bottom and the top of the submerged object that acts on its surfaces in a fluid [52, 53], and it is expressed on (13), where Vsubmerged is the buoy submerged portion volume in the sea and ρw is the water density [49, 54].(13)Fbuoyancy=ρw⋅g⋅Vsubmerged.In this case, the submerged volume change is influenced by the buoy position change. The volume of the submerged portion of the buoy is expressed as follows:(14)Vsubmerged=π⋅R2⋅Zwave−z,zwave−h≤z≤zwave,where h and R are cylindrical height and radius, respectively, Zwave is the wave height, and z coordinate represents the buoy depth from the perspective of the bottom surface [53, 54]. If z=0, the buoy bottom surface is level to ocean top surface; if z=−6, the buoy is submerged 6 meters below the ocean surface. In the case of spherical shape, the submerged volume is expressed in(15)Vsubmerged=43⋅π⋅R3⋅Zwave−z.The volume of the cylindrical buoy is determined by its radius and height, while the spherical shaped buoy volume is only determined by its radius. These geometric construction characteristics can have a different influence on the energy extraction of the WEC system. With that factor, the experiment will determine which buoy shape has the upper hand on wave energy extraction.


## 4.3. Fluid Drag Force

Drag force becomes present when the body and a fluid interact to produce the relative motion, and it acts opposite to body motion direction [49, 52]. The selected expression of the drag force is based on the Reynolds number [44], Re, which computes the flow condition, as expressed on (16):(16)Re=ρw⋅v⋅Lμw,where μw is the dynamic water viscosity, v is the relative buoy velocity from the water perspective, and L is the characteristic linear dimension [55]. Drag force is expressed as follows:(17)Fdrag=12ρw⋅v2⋅Cd⋅A,where Cd is the drag coefficient and A is the surface area influenced by the drag force [48, 50].


## 4.4. Excitation Force

The most intricate force exerted on the buoy is the wave excitation force. It is a superposition of the wave’s incident, diffraction, and radiation potentials [56]. If the excitation force is assumed to be equivalent to the incident potential, it is referred to as the Froude-Krylov force [39, 56], expressed on (18).(18)FFK=−ρw∬sub∂ΦI∂tn^⋅dS.The partial derivative from the perspective of time is done for the potential of the incident waveΦI, while the surface integral is taken over the submerged surfaces, where n^ is a normal vector unit to the surface, and dS is the dot product with the incremental surface area [49, 52]. The incident potential is expressed on (19), where A is the incident wave amplitude, ωw is the angular frequency in radian per second, k is the angular wavenumber in radian per meter, z is the depth of the submerged buoy, t is time, x is the horizontal position, and Re is the operator that keeps the actual component of a complex number [39].(19)ΦI=A⋅ωwkekzReieiωwt−kx.After the Froude-Krylov formula is simplified, the derived excitation force is estimated as in (20), where vwave is the vertical velocity of the wave [39].(20)FEXC=ρw⋅πR2k⋅vwave⋅k⋅z.


## 4.5. Mechanical System Model

The model of a mechanical system combines all the forces that act on the cable-driven pulley. Newton’s second law of motion is used to express the rotational motion, as expressed in (21), where I is the moment of inertia of the whole system being rotated by the pulley, α is the pulley angular acceleration, and τpulley are the torques that are laid on the pulley [49, 50].(21)I⋅a=∑τpulley.Equation (22) shows the resistive torque τres, where βback and βfric are the coefficients used for back-torque and mechanical friction, respectively, and CLoad is the coefficient of the load [49].(22)τres=CLoad⋅βback+βfric⋅ω.


## 4.6. Power Output and Per-Unit Cost

The WEC system output power from the linear generator is approximated to the square of the rotor’s angular velocity [49]. The instantaneous power P, generated by the generator, is given by (23), where Cpower is the generator-dependent power coefficient [49].(23)P=CLoad⋅Cpower⋅ω2.Objective one is derived based on (23), where fx1 is equated to the sum of the power generated by the liar generator. The WEC system’s annual generated energy is expressed on (24). During the multiobjective optimization of the WEC system, objective one will be maximized while objective 2 as expressed on (27) will be minimized.(24)Objective1:fx1=∑Pave,were Pave is the average power generated by a linear generator that is expected to be maximized.


## 4.7. Annual Energy

(25)Eannual=P⋅Hoursdaily⋅daysyearly.LCOE, expressed on (26), is a measure of the average net present cost of energy for a generating plant lifetime [54]. In this case, the generating period is set to 1 year to determine the average cost of the per-unit energy vs. the total annual power generated by the optimized WEC system.(26)LCOE=Eannualm−0.5,(27)Objective2:fx2=LCOE,where LCOE is the Levelized Cost of Energy to be minimized, Eannual is the total annual energy, m is the buoy mass in kilograms, and −0.5 is inserted for the optimization process based on [54].


## 5. Experimental Discussion

The benchmark functions F1, F2, F3, F4, F5, F6, and F7, which have complicated Pareto fronts, are employed to examine the proposed HMOONDS-CD algorithm’s performance [5], and it is compared to other algorithms’ performance, namely, the original MOEA/D [19], the improved iMOEA/D [5], NSGA-II [4], MOSMA [11], and MOMVO [57]. The Matlab 2017 software is employed for the experiment simulations. The proposed HMOONDS-CD algorithm is further applied to the WEC multiobjective optimization problem, and the performance is compared to other state-of-the-art multiobjective algorithms. The specifications of the hardware used are Lenovo G40, Processor Intel core (TM) i5-4210U CPU@1.70 GHz 2.40 GHz, RAM installed 4 GB, System Type 64-bit operating system, and x 64-based processor. The performance matrices indicators are employed to examine the algorithms’ performance [58, 59]. Two performance indicators with different characters are employed for a fair and unbiased comparison of the algorithms, namely, the HV and IGD. The below discussion outlines the difference between these performance indicators.


### 5.1. Hypervolume Indicator (HV)

WithZ∗=Z1∗,…,Zm∗ being the reference point in the region of the objective which fulfils Zj∗≥maxfj and P being an estimate set of PF profited by an algorithm [60, 61], the hypervolume of P is the volume of the dominated search space by P and is bounded by Z∗ and is determined by(28)HVP,Z∗=volumeUf1x,Z1∗×⋯fmx,Zm∗,where volume (.) is the Lebesgue measure, and the algorithm that achieves the largest HV metric is considered the better algorithm [5]. The hypervolume performance indicator is mainly employed for the PF’s convergence and distribution performance indicator [62].


### 5.2. Inverted Generational Distance (IGD)

For the IGD definitionP∗ is set to be the set of points uniformly distributed over the PF in the region of the objective, and like on hypervolume indicator, P is set to be an estimate set of PF profited by an algorithm [5, 61, 63, 64]. Therefore the inverted generational distance from P∗ to P is expressed as follows:(29)IGDP∗,P=∑p∈PdV,P∗P.dV,P∗ is the Euclidean distance between the member V of P and the nearest member of P∗ [5]. The algorithm that achieves the lowest IGD metric is considered the better algorithm [5].


### 5.3. Wilcoxon Signed-Rank Test

Wilcoxon signed-rank test is used to perform the nonparametric statistical analyses for this study. The Wilcoxon signed-rank test is a nonparametric test that compares the data with a population that does not have a normal distribution. A level of significance ofa = 0.05 is employed during the Wilcoxon test, where a P value indicates whether a test is significant or not, and it also reflects how significant the result is [65]. A smaller P value than a = 0.05 indicates a more vital significance, and a null hypothesis is accepted, while a larger P value than a = 0.05 indicates a weaker significance, and the null hypothesis is rejected, and an alternative hypothesis is considered.


### 5.4. Benchmark Function Problems Results

Each function will be run independently 30 times. Table1 gives the parameter settings for the algorithms to be compared on this paper.Table 1
Parameter settings of the algorithms.
AlgorithmNMaxgenTmTrmaxδFCrC1C2ωMEA/D1005000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]———iMEA/D1005000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]———HMOONDS-CD1005000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]22[0.8 1]NSGA-II100500————————MOSMA100500————————MOMVO100500————————Based on Pareto front (PF) results of IGD metric and HV metric shown on Tables2 and 3, respectively, it can be noticed for F1 problem that the improved HMOONDS-CD performed better than the NSGA-II, iMOEA/D, MOSMA, and MOMVO on IGD metrics PF, with their best IGD metric values on simulations being 0.0312, 0.0081, 0.8259, and 0.1744 respectively, which are very high compared to 0.0067 of improved HMOONDS-CD. The original MOEA/D performed slightly better than the improved HMOONDS-CD, making the best metric value of 0.0064, which is just 3% better than the improved HMOONDS-CD. However, on HV metric results for F1 problem, the improved HMOONDS-CD outperformed the other algorithms. Its best HV metric value is 1.9219, which is far better than the best values of NSGA-II, iMOEA/D, MOEA/D, MOSMA, and MOMVO, which are 1.6297, 1.7342, 1.6126, 0.2452, and 0.9953, respectively. The rest of the results are shown on Tables 2 and 3, with the best performance values in bold.Table 2
Inverted generational distance metrics.
ProblemMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOF1Best0.00640.00810.03120.00670.82590.1744Mean0.007220.009780.034990.00778330.99980.2182Worse0.00870.01120.04410.00961.27670.2818F2Best0.0010.0010.00610.0011.24660.1883Mean0.001510.001330.011030.001271.44140.2670Worse0.00130.00170.0210.00171.50010.4558F3Best0.00250.00270.01060.00110.60310.0473Mean0.0030270.003270.0133870.001330.66400.0633Worse0.00350.00410.01630.00170.71600.1185F4Best0.00110.00340.00350.0004790.00590.0325Mean0.0011500.003650.00400.000744380.03320.0471Worse0.00130.0040.00470.000966150.10960.0817F5Best0.000747150.000898950.000939140.000617240.00680.0299Mean0.000778360.000957030.00114460.000733940.04850.0540Worse0.000831110.0010.00140.000839710.14900.0955F6Best0.02150.01170.00880.000833220.00680.0217Mean0.022350.01260.014120.00397110.02180.0417Worse0.02320.01340.01850.00790.05090.0871F7Best0.000392030.000440140.000624480.000433830.00480.0695Mean0.000458120.0005220.000744620.000623330.05990.3230Worse0.000516980.000657560.000979430.00096290.21220.3561The bold values are the best values among different algorithms.Table 3
Hypervolume metrics.
ProblemMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOF1Best1.61261.73421.62971.92190.26361.2875Mean1.59881.70761.6031.78390.24520.9953Worse1.5821.69411.56881.73250.24250.5929F2Best3.94593.94773.93663.94430.29953.8194Mean3.94543.94713.93353.94310.23953.4842Worse3.94413.94663.92393.94110.20202.7627F3Best3.48913.49123.47933.49410.50093.4859Mean3.48683.48813.47593.49320.47673.4842Worse3.48373.48453.47133.49140.45202.7627F4Best3.8953.89723.8963.87360.97843.8669Mean3.89493.8973.89563.86580.97483.4842Worse3.89483.89693.89523.85640.96782.7627F5Best3.82483.82623.82523.82330.96373.7964Mean3.82473.82623.82493.82240.95403.7751Worse3.82473.82613.82453.82050.91163.7489F6Best3.3213.32113.32293.31390.86003.2981Mean3.30913.30213.32263.20090.85713.1277Worse3.2733.27413.32222.88480.85292.7090F7Best3.67873.67843.6793.67850.93333.5331Mean3.6783.67723.67883.67790.91353.0897Worse3.6763.66943.67843.67610.84793.0285The bold values are the best values among different algorithms.Looking at the whole simulation’s total overall performance where IGD and HV metrics are combined on Table4, the improved HMOONDS-CD algorithm improved overall performance compared to the original MOEA/D, NSGA-II, iMOEA/D, MOSMA, and the MOMVO algorithms. The improved HMOONDS-CD achieved 50% on the best metric performance indicator and 50% on the mean metric performance indicator.Table 4
Overall performance analysis.
Evaluation modeMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOIGD matrixBest310500Mean200500Worse106000Overall % best42.82%14.29%0%71.43%0%0%Overall % mean28.57%0%0%71.43%0%0%Overall % worse14.29%0%85.71%0%0%0%HV matrixBest032200Mean032200Worse013300Overall % best0%42.82%28.57%28.57%0%0%Overall % mean0%42.82%28.57%28.57%0%0%Overall % worse0%14.29%42.8242.820%0%OverallTotal best342700Total mean232700Total worse119300Total overall % best21.43%28.57%14.29%50%0%0%Total overall % mean14.29%21.43%14.29%50%0%0%Total overall % worse7.14%7.14%64.29%21.43%0%0%Based on Wilcoxon signed-rank nonparametric test for both IGD and HV results on benchmark function on Tables5 and 6, respectively, it can be noticed that the statistical analysis test rejected the null hypotheses of the HV results that suggested that MOEA/D, NSGA-II, and iMOEA/D algorithms performed better than the improved HMOONDS-CD. Therefore, in this case, the alternative hypothesis is taken, which suggests that MOEA/D, NSGA-II, and iMOEA/D did not perform better than the improved HMOONDS-CD besides the MOSMA and the MOMVO. By looking at IGD statistical analyses, it can be noticed that the null hypothesis that suggested that the enhanced HMOONDS-CD performed better than the original MOEA/D was rejected, and again the alternative hypothesis is taken. However, the null hypothesis was accepted, which suggested that the improved HMOONDS-CD performed better than iMOEA/D, NSGA-II, MOSMA, and MOMVO.Table 5
HV Wilcoxon signed ranked test for benchmark functions.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.69500.305NoHMOONDS-CD vs. iMOEA/D0.69500.305NoHMOONDS-CD vs. NSGA-II0.50.5NoHMOONDS-CD vs. MOSMA0.99290.0071YesHMOONDS-CD vs. MOMVO0.99290.0071YesTable 6
IDG Wilcoxon test for benchmark functions.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.88100.119NoHMOONDS-CD vs. iMOEA/D0.97880.0212YesHMOONDS-CD vs. NSGA-II0.99110.0089YesHMOONDS-CD vs. MOSMA0.99290.0071YesHMOONDS-CD vs. MOMVO0.99290.0071Yes


### 5.5. Multiobjective Wave Energy Converter Problem Results

During the simulation of the optimization of multiobjective WEC problem, four types of problems will be simulated using the improved HMOONDS-CD in comparison with other abovementioned multiobjective algorithms. The first WEC problem type will consist of cylindrical (π.R2.h) buoy at depth z=0, meaning the cylindrical buoy bottom surface will be in the level of ocean top surface. The second WEC problem type will consist of cylindrical (π.R2.h) buoy at depth z=−6, meaning the whole cylindrical buoy will be submerged to a low level of 6 meters below the ocean top surface. The third WEC problem will consist of spherical (4/3.π.R3) buoy at depth z=0, while the fourth WEC problem consists of spherical (4/3.π.R3) buoy at depth z=−6.Following are the constant parameters of the WEC problem [39]:m=500, buoy total mass (kg).g=9.81, gravitational acceleration (m/s2).I=0.04, flywheel moment of inertia (kg-m2).To=10, tension in reel (kg-m/s2).r=0.05, radius of pulley (m).βfric=0.5, coefficient for the mechanical friction.βback=0.343, velocity-dependent parameter of back-torque (kg-m/s).Cpower=4.10, ratio of power output to the square of omega (kg-m2/s).ρw=997, density of water (kg-m3).Cd=0.82, coefficient of drag of buoy.Following are varying parameters of the WEC problem:R=X1, buoy radius (m).X1 upper bound = 5.X1 lower bound = 3.h=X2, buoy cylindrical height (m).X2 upper bound = 4.X2 lower bound = 2.Each WEC problem will be run independently 30 times. Table7 gives the parameter settings for the algorithms to be compared on this simulation session.Table 7
Parameter settings for the algorithms.
AlgorithmNMaxgenTmTrmaxδFCrC1C2ωMEA/D805000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]———iMEA/D805000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]2——HMOONDS-CD805000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]22[0.8 1]NSGA-II80500————————MOSMA80500————————MOMVO80500————————The results of cylindrical buoy and spherical buoy WEC optimization problems are shown in Table8. Based on the cylindrical WEC problem results, where Z=0, it can be noticed that the improved HMOONDS-CD archived the best maximum annual power generation of 396.500 TW (Terra watts) and best minimum per-unit cost of power of 22 c (cents) compared to other algorithms. The second-best algorithm is the NSGA-II which achieved the maximum annual power generation of 340.700 TW and average power of 283.836 TW, while achieving minimum and average per-unit power cost of 32 c and 32 c, respectively. The improved HMOONDS-CD achieved an average power of 311.203 TW and an average per-unit cost of 32 c, which is best compared to the second-best algorithm of NSGA-II. The iMOEA/D achieved the most unsatisfactory results compared to all other algorithms, with the maximum annual power of 167.700 TW and minimum per unit of 34 c. MOSMA and MOMVO achieved the best maximum power of 333.860 TW and 230.341 TW, respectively, and the best-LCOE of 32 c and 33 c, respectively. Looking at the cylindrical buoy WEC problem of Z=−6, the improved HMOONDS-CD achieved the maximum and average annual generated power of 363.900 TW and 312.680 TW, respectively, the best results compared to other optimization results by other algorithms.Table 8
WEC problems performances.
ProblemMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOCylindricalZ=0Best-power186.600e + 12167.700e + 12340.700e + 12396.500e + 12333.860e + 12230.341e + 12Mean-power128.110e + 1299.196e + 12283.836e + 12311.203e + 12288.31e + 12176.189e + 12Worse-power59.2300e + 1244.430e + 12240.000e + 12243.000e + 12272.70e + 1280.280e + 12Best-LCOE0.340.340.320.220.320.33Mean-LCOE0.350.360.320.320.320.34Worse-LCOE0.380.390.330.330.330.37CylindricalZ=−6Best-power207.200e + 12216.700e + 12339.000e + 12387.200e + 12355.46e + 12209.905e + 12Mean-power130.960e + 12129.490e + 12291.603e + 12314.403e + 12308.41e + 12103.520e + 12Worse-power66.4800e + 1258.280e + 12259.500e + 12258.900e + 12258.60e + 1251.906e + 12Best-LCOE0.330.330.320.310.320.33Mean-LCOE0.350.350.320.320.320.36Worse-LCOE0.370.380.330.330.330.38SphericalZ=0Best-power199.100e + 12181.300e + 12354.200e + 12405.500e + 12318.200e + 12225.039e + 12Mean-power126.940e + 12105.773e + 12278.233e + 12306.696e + 12289.860e + 12147.758e + 12Worse-power75.430e + 1229.700e + 12226.100e + 12257.700e + 12253.420e + 1247.219e + 12Best-LCOE0.340.340.320.310.320.33Mean-LCOE0.350.360.320.320.320.35Worse-LCOE0.370.410.330.330.330.39SphericalZ=−6Best-power297.600e + 12186.400e + 12397.500e + 12405.100e + 12322.70e + 12233.908e + 12Mean-power152.479e + 12113.782e + 12293.660e + 12316.880e + 12252.40e + 12149.704e + 12Worse-power95.390e + 1220.650e + 12224.000e + 12203.600e + 12200.50e + 1239.751e + 12Best-LCOE0.320.330.310.310.320.33Mean-LCOE0.340.350.320.320.320.36Worse-LCOE0.360.370.330.330.330.39Based on the spherical buoy WEC problem atZ=0, the improved HMOONDS-CD achieved the best maximum and average annual generated power of 386.200 TW and 300.725 TW, respectively, while achieving the best minimum and average cost per unit power of 31 c and 32 c, respectively. The NSGA-II achieved the second-best results again with maximum annual generated power of 354.200 TW and a minimum per-unit cost of 32 c. The original MOEA/D and iMOEA/D recorded the maximum powers and per-unit costs of 199.100 TW and 181.300 TW, respectively, and 34 c for both algorithms. MOSMA and MOMVO achieved the best maximum power of 318.200 TW and 225.039 TW, respectively, and the best-LCOE of 32 c and 33 c, respectively. Looking at the spherical buoy WEC problem at Z = −6, the improved HMOONDS-CD recorded the best maximum and average annual generated power of 405.100 TW and 318.165 TW, respectively, the highest generated annual power compared to the optimization by other algorithms. The rest of the results are shown in Table 8 above. Table 9 shows the overall performance analysis by all the employed algorithms in this paper based on the WEC problem and Table 10 shows HV and IGD performance results for WEC problems.Table 9
Overall performance analysis.
Evaluation modeMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOAnnual powerBest000400Mean000400Worse002200Overall % best0%0%0%100%0%0%Overall % mean0%0%0%100%0%0%Overall % worse0%0%50%50%0%0%Per-unit costBest001400Mean004440Worse004440Overall % best0%0%25%100%0%0%Overall % mean0%0%100%100%100%0%Overall % worse0%0%100%100%100%0%OverallTotal best001800Total mean004840Total worse004640Total overall % best0%0%12.5%100%0%0%Total overall % mean0%0%50%100%50%0%Total overall % worse0%0%50%75%50%0%Table 10
HV and IGD performance results for WEC problems.
MOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOCylinder Z = 0HV-mean−5.953e + 12−5.045e + 124−12.521e + 12−11.118e + 12−8.820e + 12IGD-mean13.1839e + 1218.679e + 121.414243.372e + 12109.631e + 06139.385e + 06Cylinder Z = −6HV-mean−6.112e + 12−6.756e + 124−12.628e + 12−11.955e + 12−8.880e + 12IGD-mean12.840e + 1220.119e + 121.414244.349e + 12123.555e + 06170.820e + 06Sphere Z = 0HV-mean−5.805e + 12−5.383e + 124−12.352e + 12−11.194e + 12−9.016e + 12IGD-mean13.083e + 1218.787e + 121.414242.886e + 12110.04e + 06161.261e + 06Sphere Z = −6HV-mean−6.761e + 12−7.224e + 124−12.720e + 12−11.585e + 12−9.042e + 12IGD-mean15.053e + 1224.709e + 121.414244.662e + 12123.31e + 06159.207e06Based on the Wilcoxon test for both IGD and HV results on WEC problems in Tables11 and 12, respectively, the null hypotheses were accepted based on the significance level α=0.05, which suggested that the improved HMOONDS-CD performed better than other algorithms for the HV results in Table 11 and that all other algorithms performed better than the enhanced HMOONDS-CD for the IGD results in Table 12.Table 11
HV Wilcoxon test for WEC problem.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.96560.0344YesHMOONDS-CD vs. iMOEA/D0.96560.0344YesHMOONDS-CD vs. NSGA-II0.96560.0344YesHMOONDS-CD vs. MOSMA0.96560.0344YesHMOONDS-CD vs. MOMVO0.96560.0344YesTable 12
IDG Wilcoxon test for WEC problem.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.96560.0344YesHMOONDS-CD vs. iMOEA/D0.96560.0344YesHMOONDS-CD vs. NSGA-II0.96560.0344YesHMOONDS-CD vs. MOSMA0.96560.0344YesHMOONDS-CD vs. MOMVO0.96560.0344Yes


## 5.1. Hypervolume Indicator (HV)

WithZ∗=Z1∗,…,Zm∗ being the reference point in the region of the objective which fulfils Zj∗≥maxfj and P being an estimate set of PF profited by an algorithm [60, 61], the hypervolume of P is the volume of the dominated search space by P and is bounded by Z∗ and is determined by(28)HVP,Z∗=volumeUf1x,Z1∗×⋯fmx,Zm∗,where volume (.) is the Lebesgue measure, and the algorithm that achieves the largest HV metric is considered the better algorithm [5]. The hypervolume performance indicator is mainly employed for the PF’s convergence and distribution performance indicator [62].


## 5.2. Inverted Generational Distance (IGD)

For the IGD definitionP∗ is set to be the set of points uniformly distributed over the PF in the region of the objective, and like on hypervolume indicator, P is set to be an estimate set of PF profited by an algorithm [5, 61, 63, 64]. Therefore the inverted generational distance from P∗ to P is expressed as follows:(29)IGDP∗,P=∑p∈PdV,P∗P.dV,P∗ is the Euclidean distance between the member V of P and the nearest member of P∗ [5]. The algorithm that achieves the lowest IGD metric is considered the better algorithm [5].


## 5.3. Wilcoxon Signed-Rank Test

Wilcoxon signed-rank test is used to perform the nonparametric statistical analyses for this study. The Wilcoxon signed-rank test is a nonparametric test that compares the data with a population that does not have a normal distribution. A level of significance ofa = 0.05 is employed during the Wilcoxon test, where a P value indicates whether a test is significant or not, and it also reflects how significant the result is [65]. A smaller P value than a = 0.05 indicates a more vital significance, and a null hypothesis is accepted, while a larger P value than a = 0.05 indicates a weaker significance, and the null hypothesis is rejected, and an alternative hypothesis is considered.


## 5.4. Benchmark Function Problems Results

Each function will be run independently 30 times. Table1 gives the parameter settings for the algorithms to be compared on this paper.Table 1
Parameter settings of the algorithms.
AlgorithmNMaxgenTmTrmaxδFCrC1C2ωMEA/D1005000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]———iMEA/D1005000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]———HMOONDS-CD1005000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]22[0.8 1]NSGA-II100500————————MOSMA100500————————MOMVO100500————————Based on Pareto front (PF) results of IGD metric and HV metric shown on Tables2 and 3, respectively, it can be noticed for F1 problem that the improved HMOONDS-CD performed better than the NSGA-II, iMOEA/D, MOSMA, and MOMVO on IGD metrics PF, with their best IGD metric values on simulations being 0.0312, 0.0081, 0.8259, and 0.1744 respectively, which are very high compared to 0.0067 of improved HMOONDS-CD. The original MOEA/D performed slightly better than the improved HMOONDS-CD, making the best metric value of 0.0064, which is just 3% better than the improved HMOONDS-CD. However, on HV metric results for F1 problem, the improved HMOONDS-CD outperformed the other algorithms. Its best HV metric value is 1.9219, which is far better than the best values of NSGA-II, iMOEA/D, MOEA/D, MOSMA, and MOMVO, which are 1.6297, 1.7342, 1.6126, 0.2452, and 0.9953, respectively. The rest of the results are shown on Tables 2 and 3, with the best performance values in bold.Table 2
Inverted generational distance metrics.
ProblemMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOF1Best0.00640.00810.03120.00670.82590.1744Mean0.007220.009780.034990.00778330.99980.2182Worse0.00870.01120.04410.00961.27670.2818F2Best0.0010.0010.00610.0011.24660.1883Mean0.001510.001330.011030.001271.44140.2670Worse0.00130.00170.0210.00171.50010.4558F3Best0.00250.00270.01060.00110.60310.0473Mean0.0030270.003270.0133870.001330.66400.0633Worse0.00350.00410.01630.00170.71600.1185F4Best0.00110.00340.00350.0004790.00590.0325Mean0.0011500.003650.00400.000744380.03320.0471Worse0.00130.0040.00470.000966150.10960.0817F5Best0.000747150.000898950.000939140.000617240.00680.0299Mean0.000778360.000957030.00114460.000733940.04850.0540Worse0.000831110.0010.00140.000839710.14900.0955F6Best0.02150.01170.00880.000833220.00680.0217Mean0.022350.01260.014120.00397110.02180.0417Worse0.02320.01340.01850.00790.05090.0871F7Best0.000392030.000440140.000624480.000433830.00480.0695Mean0.000458120.0005220.000744620.000623330.05990.3230Worse0.000516980.000657560.000979430.00096290.21220.3561The bold values are the best values among different algorithms.Table 3
Hypervolume metrics.
ProblemMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOF1Best1.61261.73421.62971.92190.26361.2875Mean1.59881.70761.6031.78390.24520.9953Worse1.5821.69411.56881.73250.24250.5929F2Best3.94593.94773.93663.94430.29953.8194Mean3.94543.94713.93353.94310.23953.4842Worse3.94413.94663.92393.94110.20202.7627F3Best3.48913.49123.47933.49410.50093.4859Mean3.48683.48813.47593.49320.47673.4842Worse3.48373.48453.47133.49140.45202.7627F4Best3.8953.89723.8963.87360.97843.8669Mean3.89493.8973.89563.86580.97483.4842Worse3.89483.89693.89523.85640.96782.7627F5Best3.82483.82623.82523.82330.96373.7964Mean3.82473.82623.82493.82240.95403.7751Worse3.82473.82613.82453.82050.91163.7489F6Best3.3213.32113.32293.31390.86003.2981Mean3.30913.30213.32263.20090.85713.1277Worse3.2733.27413.32222.88480.85292.7090F7Best3.67873.67843.6793.67850.93333.5331Mean3.6783.67723.67883.67790.91353.0897Worse3.6763.66943.67843.67610.84793.0285The bold values are the best values among different algorithms.Looking at the whole simulation’s total overall performance where IGD and HV metrics are combined on Table4, the improved HMOONDS-CD algorithm improved overall performance compared to the original MOEA/D, NSGA-II, iMOEA/D, MOSMA, and the MOMVO algorithms. The improved HMOONDS-CD achieved 50% on the best metric performance indicator and 50% on the mean metric performance indicator.Table 4
Overall performance analysis.
Evaluation modeMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOIGD matrixBest310500Mean200500Worse106000Overall % best42.82%14.29%0%71.43%0%0%Overall % mean28.57%0%0%71.43%0%0%Overall % worse14.29%0%85.71%0%0%0%HV matrixBest032200Mean032200Worse013300Overall % best0%42.82%28.57%28.57%0%0%Overall % mean0%42.82%28.57%28.57%0%0%Overall % worse0%14.29%42.8242.820%0%OverallTotal best342700Total mean232700Total worse119300Total overall % best21.43%28.57%14.29%50%0%0%Total overall % mean14.29%21.43%14.29%50%0%0%Total overall % worse7.14%7.14%64.29%21.43%0%0%Based on Wilcoxon signed-rank nonparametric test for both IGD and HV results on benchmark function on Tables5 and 6, respectively, it can be noticed that the statistical analysis test rejected the null hypotheses of the HV results that suggested that MOEA/D, NSGA-II, and iMOEA/D algorithms performed better than the improved HMOONDS-CD. Therefore, in this case, the alternative hypothesis is taken, which suggests that MOEA/D, NSGA-II, and iMOEA/D did not perform better than the improved HMOONDS-CD besides the MOSMA and the MOMVO. By looking at IGD statistical analyses, it can be noticed that the null hypothesis that suggested that the enhanced HMOONDS-CD performed better than the original MOEA/D was rejected, and again the alternative hypothesis is taken. However, the null hypothesis was accepted, which suggested that the improved HMOONDS-CD performed better than iMOEA/D, NSGA-II, MOSMA, and MOMVO.Table 5
HV Wilcoxon signed ranked test for benchmark functions.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.69500.305NoHMOONDS-CD vs. iMOEA/D0.69500.305NoHMOONDS-CD vs. NSGA-II0.50.5NoHMOONDS-CD vs. MOSMA0.99290.0071YesHMOONDS-CD vs. MOMVO0.99290.0071YesTable 6
IDG Wilcoxon test for benchmark functions.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.88100.119NoHMOONDS-CD vs. iMOEA/D0.97880.0212YesHMOONDS-CD vs. NSGA-II0.99110.0089YesHMOONDS-CD vs. MOSMA0.99290.0071YesHMOONDS-CD vs. MOMVO0.99290.0071Yes


## 5.5. Multiobjective Wave Energy Converter Problem Results

During the simulation of the optimization of multiobjective WEC problem, four types of problems will be simulated using the improved HMOONDS-CD in comparison with other abovementioned multiobjective algorithms. The first WEC problem type will consist of cylindrical (π.R2.h) buoy at depth z=0, meaning the cylindrical buoy bottom surface will be in the level of ocean top surface. The second WEC problem type will consist of cylindrical (π.R2.h) buoy at depth z=−6, meaning the whole cylindrical buoy will be submerged to a low level of 6 meters below the ocean top surface. The third WEC problem will consist of spherical (4/3.π.R3) buoy at depth z=0, while the fourth WEC problem consists of spherical (4/3.π.R3) buoy at depth z=−6.Following are the constant parameters of the WEC problem [39]:m=500, buoy total mass (kg).g=9.81, gravitational acceleration (m/s2).I=0.04, flywheel moment of inertia (kg-m2).To=10, tension in reel (kg-m/s2).r=0.05, radius of pulley (m).βfric=0.5, coefficient for the mechanical friction.βback=0.343, velocity-dependent parameter of back-torque (kg-m/s).Cpower=4.10, ratio of power output to the square of omega (kg-m2/s).ρw=997, density of water (kg-m3).Cd=0.82, coefficient of drag of buoy.Following are varying parameters of the WEC problem:R=X1, buoy radius (m).X1 upper bound = 5.X1 lower bound = 3.h=X2, buoy cylindrical height (m).X2 upper bound = 4.X2 lower bound = 2.Each WEC problem will be run independently 30 times. Table7 gives the parameter settings for the algorithms to be compared on this simulation session.Table 7
Parameter settings for the algorithms.
AlgorithmNMaxgenTmTrmaxδFCrC1C2ωMEA/D805000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]———iMEA/D805000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]2——HMOONDS-CD805000.2 N0.2 N0.9[0.4 0.6][0.4 0.6]22[0.8 1]NSGA-II80500————————MOSMA80500————————MOMVO80500————————The results of cylindrical buoy and spherical buoy WEC optimization problems are shown in Table8. Based on the cylindrical WEC problem results, where Z=0, it can be noticed that the improved HMOONDS-CD archived the best maximum annual power generation of 396.500 TW (Terra watts) and best minimum per-unit cost of power of 22 c (cents) compared to other algorithms. The second-best algorithm is the NSGA-II which achieved the maximum annual power generation of 340.700 TW and average power of 283.836 TW, while achieving minimum and average per-unit power cost of 32 c and 32 c, respectively. The improved HMOONDS-CD achieved an average power of 311.203 TW and an average per-unit cost of 32 c, which is best compared to the second-best algorithm of NSGA-II. The iMOEA/D achieved the most unsatisfactory results compared to all other algorithms, with the maximum annual power of 167.700 TW and minimum per unit of 34 c. MOSMA and MOMVO achieved the best maximum power of 333.860 TW and 230.341 TW, respectively, and the best-LCOE of 32 c and 33 c, respectively. Looking at the cylindrical buoy WEC problem of Z=−6, the improved HMOONDS-CD achieved the maximum and average annual generated power of 363.900 TW and 312.680 TW, respectively, the best results compared to other optimization results by other algorithms.Table 8
WEC problems performances.
ProblemMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOCylindricalZ=0Best-power186.600e + 12167.700e + 12340.700e + 12396.500e + 12333.860e + 12230.341e + 12Mean-power128.110e + 1299.196e + 12283.836e + 12311.203e + 12288.31e + 12176.189e + 12Worse-power59.2300e + 1244.430e + 12240.000e + 12243.000e + 12272.70e + 1280.280e + 12Best-LCOE0.340.340.320.220.320.33Mean-LCOE0.350.360.320.320.320.34Worse-LCOE0.380.390.330.330.330.37CylindricalZ=−6Best-power207.200e + 12216.700e + 12339.000e + 12387.200e + 12355.46e + 12209.905e + 12Mean-power130.960e + 12129.490e + 12291.603e + 12314.403e + 12308.41e + 12103.520e + 12Worse-power66.4800e + 1258.280e + 12259.500e + 12258.900e + 12258.60e + 1251.906e + 12Best-LCOE0.330.330.320.310.320.33Mean-LCOE0.350.350.320.320.320.36Worse-LCOE0.370.380.330.330.330.38SphericalZ=0Best-power199.100e + 12181.300e + 12354.200e + 12405.500e + 12318.200e + 12225.039e + 12Mean-power126.940e + 12105.773e + 12278.233e + 12306.696e + 12289.860e + 12147.758e + 12Worse-power75.430e + 1229.700e + 12226.100e + 12257.700e + 12253.420e + 1247.219e + 12Best-LCOE0.340.340.320.310.320.33Mean-LCOE0.350.360.320.320.320.35Worse-LCOE0.370.410.330.330.330.39SphericalZ=−6Best-power297.600e + 12186.400e + 12397.500e + 12405.100e + 12322.70e + 12233.908e + 12Mean-power152.479e + 12113.782e + 12293.660e + 12316.880e + 12252.40e + 12149.704e + 12Worse-power95.390e + 1220.650e + 12224.000e + 12203.600e + 12200.50e + 1239.751e + 12Best-LCOE0.320.330.310.310.320.33Mean-LCOE0.340.350.320.320.320.36Worse-LCOE0.360.370.330.330.330.39Based on the spherical buoy WEC problem atZ=0, the improved HMOONDS-CD achieved the best maximum and average annual generated power of 386.200 TW and 300.725 TW, respectively, while achieving the best minimum and average cost per unit power of 31 c and 32 c, respectively. The NSGA-II achieved the second-best results again with maximum annual generated power of 354.200 TW and a minimum per-unit cost of 32 c. The original MOEA/D and iMOEA/D recorded the maximum powers and per-unit costs of 199.100 TW and 181.300 TW, respectively, and 34 c for both algorithms. MOSMA and MOMVO achieved the best maximum power of 318.200 TW and 225.039 TW, respectively, and the best-LCOE of 32 c and 33 c, respectively. Looking at the spherical buoy WEC problem at Z = −6, the improved HMOONDS-CD recorded the best maximum and average annual generated power of 405.100 TW and 318.165 TW, respectively, the highest generated annual power compared to the optimization by other algorithms. The rest of the results are shown in Table 8 above. Table 9 shows the overall performance analysis by all the employed algorithms in this paper based on the WEC problem and Table 10 shows HV and IGD performance results for WEC problems.Table 9
Overall performance analysis.
Evaluation modeMOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOAnnual powerBest000400Mean000400Worse002200Overall % best0%0%0%100%0%0%Overall % mean0%0%0%100%0%0%Overall % worse0%0%50%50%0%0%Per-unit costBest001400Mean004440Worse004440Overall % best0%0%25%100%0%0%Overall % mean0%0%100%100%100%0%Overall % worse0%0%100%100%100%0%OverallTotal best001800Total mean004840Total worse004640Total overall % best0%0%12.5%100%0%0%Total overall % mean0%0%50%100%50%0%Total overall % worse0%0%50%75%50%0%Table 10
HV and IGD performance results for WEC problems.
MOEA/DiMOEA/DNSGA-IIHMOONDS-CDMOSMAMOMVOCylinder Z = 0HV-mean−5.953e + 12−5.045e + 124−12.521e + 12−11.118e + 12−8.820e + 12IGD-mean13.1839e + 1218.679e + 121.414243.372e + 12109.631e + 06139.385e + 06Cylinder Z = −6HV-mean−6.112e + 12−6.756e + 124−12.628e + 12−11.955e + 12−8.880e + 12IGD-mean12.840e + 1220.119e + 121.414244.349e + 12123.555e + 06170.820e + 06Sphere Z = 0HV-mean−5.805e + 12−5.383e + 124−12.352e + 12−11.194e + 12−9.016e + 12IGD-mean13.083e + 1218.787e + 121.414242.886e + 12110.04e + 06161.261e + 06Sphere Z = −6HV-mean−6.761e + 12−7.224e + 124−12.720e + 12−11.585e + 12−9.042e + 12IGD-mean15.053e + 1224.709e + 121.414244.662e + 12123.31e + 06159.207e06Based on the Wilcoxon test for both IGD and HV results on WEC problems in Tables11 and 12, respectively, the null hypotheses were accepted based on the significance level α=0.05, which suggested that the improved HMOONDS-CD performed better than other algorithms for the HV results in Table 11 and that all other algorithms performed better than the enhanced HMOONDS-CD for the IGD results in Table 12.Table 11
HV Wilcoxon test for WEC problem.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.96560.0344YesHMOONDS-CD vs. iMOEA/D0.96560.0344YesHMOONDS-CD vs. NSGA-II0.96560.0344YesHMOONDS-CD vs. MOSMA0.96560.0344YesHMOONDS-CD vs. MOMVO0.96560.0344YesTable 12
IDG Wilcoxon test for WEC problem.
ProblemZ-scoreP valueα=0.05HMOONDS-CD vs. MOEA/D0.96560.0344YesHMOONDS-CD vs. iMOEA/D0.96560.0344YesHMOONDS-CD vs. NSGA-II0.96560.0344YesHMOONDS-CD vs. MOSMA0.96560.0344YesHMOONDS-CD vs. MOMVO0.96560.0344Yes


## 6. Conclusion

Based on the research of the improved HMOONDS-CD, it has been noticed that the algorithm was able to achieve satisfactory results during the optimization of the multiobjective problems with complex Pareto front shape. The proposed HMOONDS-CD achieved 50% of the best optimization and 50% of the best mean optimization based on the used performance metrics indicators, which is high by 21.43% for the best and by 28.57% for the best mean, compared to the second-best iMOEA/D algorithm which achieved 28.57% for the best optimization and 21.43% for the overall best mean. Based on the optimization of WEC problems, it can be noticed that the improved HMOONDS-CD outperformed other algorithms by achieving 100% best performance for the maximum and average annual power generation. The improved HMOONDS-CD achieved 100% for the minimum and average per-unit cost of power compared to other algorithms. The HMOONDS-CD also achieved 100% for the overall performance during optimization of the WEC problem. Based on Wilcoxon signed-rank test analyses, the statistical performance suggested that the improved HMOONDS-CD performed better on the HV metric for WEC problems. It was not better than other algorithms on the IGD metric performance. Therefore it can be concluded that the improved HMOONDS-CD performed better during optimization of multiobjective with characteristics of complex Pareto front, and it performed well during application of real-world problem of optimization of point absorber wave energy converter with the objective to maximize the annual generated power and objective to minimize the per-unit cost of the power. Figures4(a)–4(d) and 5(a)–5(d) show the convergence and distribution of the PF on benchmark problems and WEC optimization problems, respectively.Figure 4
(a) MOEA/D F1 function simulation, (b) HMOONDS-CD F1 function simulation, (c) NSGA-II F1 function simulation, (d) iMOEA/D F1 function simulation.
(a)(b)(c)(d)Figure 5
(a): MOEA/D WEC cylinder atZ = 0 simulation, (b) iMOEA/D WEC cylinder at Z = 0 simulation, (c) NSGA-II WEC cylinder at Z = 0 simulation, (d) HMOONDS-CD WEC cylinder at Z = 0 simulation.
(a)(b)(c)(d)


## 7. Future Works

For future works, more Multiobjective Evolutionary Algorithms can be improved for better PF convergence and distribution by using reinforcement machine learning algorithms such Q learning, Monte Carlo, SARSA, and Deep Q Learning. Machine learning algorithms can play a pivotal role in MOEAs by bringing the independent control of MOEAs control variables that can influence the MOEAs performance.



---

*Source: 8309697-2022-03-26.xml*