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

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. +e 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. +is 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. +e 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. +is 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. +e experimental results show that the proposedmethod 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.


Introduction
e 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 populationbased 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]. e 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].
ere 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]. e 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.
ere 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][13][14][15][16][17][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]. e 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]. e 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. ese stats indicate how much of a research magnet the WEC topic is to the industries and research institutions. is 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. erefore, some means of power compensating systems are needed in place to achieve a smooth power output [30]. ese 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. is 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/m 2 compared to other renewable sources such as solar and wind energy, which both have a low average power density of 1 kW/m 2 each, at peak solar insolation and 12 m/s, respectively [32].
is 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. e point absorber benefits are the low complexity and high reliability, which are significant for an offshore platform, as stated by [30]. is paper proposes to optimize the point absorber WEC device with the objectives of maximizing the annual generated power and minimizing the cost of the perunit 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. e 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]. eir 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]. eir study investigated a novel design of the WEC harvester that utilized a 2 International Transactions on Electrical Energy Systems 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]. e 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. e MOPs with complicated Pareto optimal front will be used to evaluate the performance of the proposed improved algorithm. e improved algorithm will further be applied to the point absorber WEC optimization problem. e 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. e section also discusses the framework of the improved MOEA/D with its fixture. Section 4 discusses the experiment setup and the performance matrix indicators to be used in this research. e 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.

Multiobjective Optimization Problem and
Original MOEA/D Background

Multiobjective Optimization Problem.
A minimized multiobjective optimization problem (MOP) is defined as follows: where x � x 1 , x 2 , . . . , x n ∈ Ω 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. R m 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]. is 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: Let q, p ∈ Ω, be the two solutions, q is said to dominate p if, and only if, q i ≤ p i for every i ∈ 1, . . . , m { } and q j < p j for at least one index j ∈ 1, . . . , m { }, denoted as q < p [1,38]. Condition 2: A feasible solution x * ∈ Ω of problem (1) is said to be Pareto optimal if there is no point x ∈ Ω such that F(x) dominates F(x * ).F(x * ) and it is called Pareto optimal (objective) vector [19]. Condition 3: e set of all the Pareto optimal points is called the Pareto optimal set P * , denoted as . Condition 4: e P * image in the objective region is called the Pareto front PF * , denoted as

Multiobjective Evolutionary Algorithm Based on Decomposition Framework.
e 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]. e 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]. e preferred information is represented based on user-specified inputs [40]. One of the preferred information is the ideal reference solution (Z ref ) to the MOP [40]. e other preferred information is the weight vector (λ i ) that assigns the relative preference of the user [40]. e objective function of the subproblem to be minimized is defined as follows [19]: where g te 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. e 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 B(i) [19]. f j (x) is the jth objective [20].
e 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 International Transactions on Electrical Energy Systems λ 1 , . . . , λ N is a group of even distributed weight vectors and Z * is the reference solution. e problem on (1) becomes the PF approximation problem that is decomposed into N of single-objective optimization subproblems. e j th subproblem's objective function is given as follows: where In this case, all the N objective functions are optimized simultaneously in a single run, and λ is influential to g te . At this instant, with the scalar objective function g te (x|λ i , Z * ) the optimal solution is relatively close to the optimal solution of the scalar objective function g te (x|λ j , Z * ) when their weight vectors λ i and λ j are close to each other. erefore, other scalar objective functions (g te (x|λ 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.
e i th subproblem's neighborhood is made up of all the subproblems with the weight vectors from the neighborhood ofλ i . e 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].

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. erefore 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. e 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] Z nad , and the initialization of the subproblems.

e Nondominated Sorting and Crowding Distance Application on the Proposed HMOONDS-CD.
In the improved MOEA/D, for each solution on the population P space Ω and the Fitness Value FV space Ω , the NdS will compute the number of all the solutions N p which dominate the solution P and the number of all the solutions N FV that dominate the solution FV. NdS will also compute the groups of solutions G p and G FV that the solution P and solution FV dominate, respectively [2,42]. e 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). erefore, for each solution P with N p � 0 and each solution FV with N FV � 0, their members q p and q FV of their groups G p and G FV 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. e same procedure is continued for identifying the third member and until all the members are identified. e CD will approximate the density of the solutions from the NdS in P ∈ Ω and FV ∈ Ω. is 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].

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. is 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. e 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. e 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]. e 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 r 1 , r 2 , r 3 (r 1 ≠ r 2 ≠ r 3 ≠ i) are randomly selected from a set of subproblems B i in order to generate the new solution y. DE/ Rand/1 mutation strategy is selected, and it is expressed as follows: where r 1 , r 2 , r 3 are three individuals randomly chosen from the set B i while X r2 ,-X r3 , are the different vectors to mutate the parent and X r1 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: 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]. e new solution obtained from DE mutation is further updated based on the Swarm Velocity and position update on the PSO phase. is strategy enhances the algorithm optimization speed for fast selection of the best possible solutions to be fed to the decomposition stage. e 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. e new updated y from PSO operation is established based on the velocity swarm update where c 1 and c 2 are learning factor 1 and global learning factor 2, respectively [47]. Parameters r 1 and r 2 are random parameters with r 1 ranging between [pBest, y] values while r 2 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. e 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. index

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]. e high cost of wave energy limits it to International Transactions on Electrical Energy Systems 5 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. is 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. e 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].
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 Figure 3 [49,50].

Buoy Weight.
e weight being the simplest force that acts on the buoy, it is expressed on (12), where m buoy is the mass of the entire buoy system, and g is the gravitational constant [49,51].
e 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].

Buoyancy Force.
e buoyancy force is exerted by the fluid in the upward direction of a partially or fully submerged object, opposing the object's weight. e 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 V submerged is the buoy submerged portion volume in the sea and ρ w is the water density [49,54].
In this case, the submerged volume change is influenced by the buoy position change. e volume of the submerged portion of the buoy is expressed as follows: where h and R are cylindrical height and radius, respectively, Z wave 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 e volume of the cylindrical buoy is determined by its radius and height, while the spherical shaped buoy volume is only determined by its radius. ese 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.

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]. e selected expression of the drag force is based on the Reynolds number [44], Re, which computes the flow condition, as expressed on (16): 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: where C d is the drag coefficient and A is the surface area influenced by the drag force [48,50].

Excitation
Force. e 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).
e 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]. e 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].
After the Froude-Krylov formula is simplified, the derived excitation force is estimated as in (20), where v wave is the vertical velocity of the wave [39].

Mechanical System Model.
e 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].
Equation (22) shows the resistive torque τ res , where β back and β fric are the coefficients used for back-torque and mechanical friction, respectively, and C Load is the coefficient of the load [49].

Power Output and Per-Unit
Cost. e WEC system output power from the linear generator is approximated to the square of the rotor's angular velocity [49]. e instantaneous power P, generated by the generator, is given by (23), where C power is the generator-dependent power coefficient [49].
Objective one is derived based on (23), where f 1 (x) is equated to the sum of the power generated by the liar generator. e 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.
were P ave is the average power generated by a linear generator that is expected to be maximized.

Annual Energy.
E annual � P · Hours daily · days yearly .
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.
Objective 2: where LCOE is the Levelized Cost of Energy to be minimized, E annual is the total annual energy, m is the buoy mass in kilograms, and −0.5 is inserted for the optimization process based on [54].

Experimental Discussion
e 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]. e Matlab 2017 software is employed for the experiment simulations. e 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. e 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. e 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 International Transactions on Electrical Energy Systems the algorithms, namely, the HV and IGD. e below discussion outlines the difference between these performance indicators.

Hypervolume Indicator (HV).
With Z * � (Z * 1 , . . . , Z * m ) being the reference point in the region of the objective which fulfils Z * j ≥ max(f j ) 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 where volume (.) is the Lebesgue measure, and the algorithm that achieves the largest HV metric is considered the better algorithm [5]. e hypervolume performance indicator is mainly employed for the PF's convergence and distribution performance indicator [62].

Inverted Generational Distance (IGD).
For the IGD definition P * 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]. erefore the inverted generational distance from P * to P is expressed as follows: d(V, P * ) is the Euclidean distance between the member V of P and the nearest member of P * [5]. e algorithm that achieves the lowest IGD metric is considered the better algorithm [5].

Wilcoxon Signed-Rank Test.
Wilcoxon signed-rank test is used to perform the nonparametric statistical analyses for this study. e 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 of a � 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.

Benchmark Function Problems Results
. Each function will be run independently 30 times. Table 1 gives the parameter settings for the algorithms to be compared on this paper. Based on Pareto front (PF) results of IGD metric and HV metric shown on Tables 2 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. e 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. e rest of the results are shown on Tables 2 and 3, with the best performance values in bold.
Looking at the whole simulation's total overall performance where IGD and HV metrics are combined on Table 4, the improved HMOONDS-CD algorithm improved overall performance compared to the original MOEA/D, NSGA-II, iMOEA/D, MOSMA, and the MOMVO algorithms. e improved HMOONDS-CD achieved 50% on the best metric performance indicator and 50% on the mean metric performance indicator.
Based on Wilcoxon signed-rank nonparametric test for both IGD and HV results on benchmark function on Tables 5 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. erefore, 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.

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. e first WEC problem type will consist of cylindrical (π.R 2 .h) buoy at depth z � 0, meaning the cylindrical buoy bottom surface will be in the level of ocean top surface. e second WEC problem type will consist of cylindrical (π.R 2 .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. e third WEC problem will consist of spherical ((4/3).π.R 3 ) buoy at depth z � 0, while the fourth WEC problem consists of spherical ((4/3).π.R 3 ) buoy at depth z � −6.
Following are varying parameters of the WEC problem: R � X 1 , buoy radius (m). X 1 upper bound � 5.
Each WEC problem will be run independently 30 times. Table 7 gives the parameter settings for the algorithms to be compared on this simulation session. e results of cylindrical buoy and spherical buoy WEC optimization problems are shown in Table 8. 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. e second-best algorithm is the NSGA-II which achieved the maximum annual power generation of 340.700 TW and average power of  -II  100  500  --------MOSMA  100  500  --------MOMVO  100  500  -------- International Transactions on Electrical Energy Systems 9 283.836 TW, while achieving minimum and average per-unit power cost of 32 c and 32 c, respectively. e 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. e iMOEA/D achieved the most unsatisfactory results compared to all other      Based on the spherical buoy WEC problem at Z � 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. e 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. e 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. e 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.      Based on the Wilcoxon test for both IGD and HV results on WEC problems in Tables 11 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.

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. e 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. e improved HMOONDS-CD achieved 100% for the minimum and average per-unit cost of power compared to other algorithms. e 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. erefore 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.

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.

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.