Constraint handling techniques for metaheuristics: a state-of-the-art review and new variants

Metaheuristic optimization algorithms (MOAs) are computational randomized search processes which draw inspiration from physical and biological phenomena, with an application spectrum that extends to numerous fields, ranging from engineering design to economics. MOAs were originally developed for solving unconstrained NP-complete problems, and hence their application to constrained optimization problems (COPs) requires the implementation of specialized techniques that facilitate the treatment of performance and bound constraints. While considerable research efforts have been oriented towards the development and subsequent enhancement of novel constraint handling techniques (CHTs) for MOAs, a systematic review of such techniques has not been conducted hitherto. This work presents a state-of-the-art review on CHTs used with MOAs and proposes eight novel variants based on the feasibility rules and ε-constrained techniques. The distinctive feature of the new variants is that they consider the level and number of constraint violations, besides the objective function value, for selection of individuals within a population. The novel variant performance is evaluated and compared with that of four well-known CHTs from the literature using the metaheuristic pity beetle algorithm, based upon 20 single-objective benchmark COPs. The computational results highlight the accuracy, effectiveness, and versatility of the novel variants, as well as their performance superiority in comparison with existing techniques, stemming from their distinctive formulation. The complete code can be downloaded from GitHub (https://github.com/nikoslagaros/MOAs-and-CHTs).


Introduction
Most modern optimization problems are pronouncedly nonlinear and contain variables subject to several function and bound constraints.Such problems are called constrained optimization problems (COPs) or nonlinear programming (NLP) problems and their formulation can be described in the following general form, according to Deb (2000): Minimi ze f (x), Subject : where f (x) is the objective function of the vector variable x, g j (x) are inequality and h k (x) are equality constraints, both referred to also as performance constraints, and x l i , x u i are respectively the lower and upper limit values of component x i in x.Equality constraints can be converted into inequality constraints as g k (x) = |h k (x)| − e ≤ 0, k = 1, 2, ..., K , where e is a small positive quantity (e.g.1.0E-3), and hence the two groups of performance constraints in expression (1) can be unified into a single group g j (x) ≤ 0, j = 1, 2, ..., m, where m = J + K .
Metaheuristic optimization algorithms (MOAs) are computational randomized search processes which draw inspiration from physical and biological phenomena.Over the years, the application field of MOAs has expanded vastly beyond physical and biological sciences and into different fields, ranging from engineering to economics.Nonetheless, most MOAs were originally developed for solving unconstrained NPcomplete optimization problems, and hence their application to COPs is rendered a challenging undertaking due to the posed constraint handling requirement.However, most of the problems where MOAs are used concern COPs.The effective overcoming of the constraint handling obstacle has progressively evolved into an independent research field, particularly since a unique constraint handling technique (CHT) which guarantees search algorithm convergence to the global optimum does not exist.The focus of research efforts has been shared between development of novel CHTs and performance enhancement of existing CHTs.With focus on MOAs, several authors have developed, enhanced, and evaluated different CHTs implemented in specific MOAs, for application to COPs across a broad range of fields.Coello and Montes (2002) conducted a detailed review and performance evaluation of the most popular CHTs implemented in evolutionary algorithms (EAs), including penalty, separation of constraints and objectives, special operators, hybrid and repair-algorithms-based CHTs.The authors further evaluated several penalty-based techniques implemented in a genetic algorithm (GA).A detailed review on several CHTs integrated into nature-inspired algorithms, e.g.EAs, or based on swarm intelligence, was conducted by Mezura-Montes and Coello (2011), including feasibility rules, ε-constrained, penalty, stochastic ranking, special operators, multi-objective concepts and ensemble CHTs.Focusing on the solution of optimal reactive power dispatch (ORPD) problems, Mallipeddi et al. (2012) evaluated penalty and adaptive penalty, superiority of feasible solutions, ε-constraint method, stochastic ranking, and ensemble CHTs, implemented in differential evolution (DE) algorithms.With focus on particle swarm optimization (PSO), the applicability and performance of penalty-based CHTs, separatists, hybrids of PSO algorithm with other optimization techniques, as well as other methods not falling under the aforementioned categories was assessed by Jordehi (2015).Miranda-Varela and Mezura-Montes (2018) evaluated four surrogate-assisted differential evolution with combined variants (SA-DECV) algorithms combined with feasibility rules, ε-constrained method, stochastic ranking, and diversity maintenance CHTs, on the basis of 24 well known test functions, and conducted comparisons with other surrogate-assisted algorithms.Ameca-Alducin et al. ( 2018) evaluated a DE algorithm combined with penalty function, feasibility rules, ε-constrained method and stochastic ranking CHTs, based upon benchmark dynamic COPs.Lin et al. (2019) assessed the performance of four different dynamic multimodal population-based optimization algorithms with superiority of feasible solution, ε-constraint method, penalty function, dynamic penalty function and stochastic selection and ranking CHTs, in the context of dynamic multimodal COPs.Caraffini et al. (2019) analyzed the behavior of various popular DE schemes combined with different CHTs, such as penalty function and correction techniques (saturation, toroidal), with focus on establishing the mutation and crossover operators that introduce structural bias in DE algorithms and appropriate combinations of CHTs, DE control parameters and population size to moderate the bias.
Notwithstanding the substantial research outlined above, a systematic review of existing CHTs for MOA has not been conducted hitherto.Within this backdrop, this paper aims to provide a detailed and complete state-of-the-art review of the most widely employed and efficient existing CHTs for MOAs.Moreover, the paper proposes eight novel variants of the well-known feasibility rules and ε-constrained CHTs and assesses their performance in comparison with existing CHTs using a swarm intelligence pity beetle algorithm (PBA) (Kallioras et al. 2018), based upon 14 mathematical and 6 engineering COPs.All computations are performed using the MATLAB variant of High-Performance Optimization Computing Platform (HP-OCP), the evolution of the original OCP (Lagaros 2014); OCP is founded on an object-oriented general-purpose code in C#, specifically developed for structural design optimization in civil engineering applications.
The paper proceeds with providing a comprehensive state-of-the-art review of almost 60 CHT studies from the literature in Sect. 2. In Sect.3, an overview of the penalty methods, feasibility rules, ε-constrained, stochastic ranking, and ensemble CHTs is provided and eight novel variants of the feasibility rules and ε-constrained CHTs are introduced.Subsequently, a brief description of HP-OCP and PBA is provided in Sect. 4. Finally, a comparative performance evaluation of the existing and novel CHTs presented in Sect. 3 is conducted in Sect.5, based upon 14 mathematical and 6 engineering COP numerical examples, followed by concluding remarks in Sect.6.

Comprehensive review on constraint handling techniques
This section provides a comprehensive review on the most widely employed and efficient existing CHTs for MOAs, classified into six categories: (i) penalty-based CHTs; (ii) CHTs based on separation of objective and constraints; (iii) combination of CHTs; (iv) repair-algorithm-based CHTs; (v) boundary-based CHTs; and (vi) other CHTs.Table 1 provides a summary of all reviewed CHTs, along with the MOA employed in the respective work and the corresponding literature reference number.

Penalty-based constraint handling techniques
Objective function penalization is the oldest and yet most widely employed CHT.Several penalty-based techniques have been developed over the years, including, among others, the static, dynamic, adaptive and death penalty methods.
During the early 2000s, Miettinen et al. (2003) tested five penalty-based CHTs with genetic algorithms to solve 33 mathematical test problems of different types, concluding that the adaptive penalty method (APM) was the most efficient and the parameter free penalty method the most reliable among them.
Nearly a decade later, Da Silva et al. ( 2011) proposed a parameter-free APM for use with DUVDE algorithm, a DE that incorporates a mechanism named Dynamic Use of Variants (DUV) that tries and automatically selects the variant with the best performance during the procedure, in the context of structural and mechanical engineering applications.This CHT, by using the feedback received from the current status of candidate solution population to adaptively define the corresponding penalty factor for every constraint combined with DUV, made DE more efficient in solving COPs and achieved very competitive results compared with GA and simple DE.Montemurro et al. (2013) proposed a parameter-free penalty-based CHT called automatic dynamic penalization (ADP) and applied it to practical engineering design problems, e.g.maximization of the first buckling load of composite laminates.The proposed CHT is combined with a multi-population genetic algorithm named BIANCA, where simultaneous evolution of different species individuals takes place by utilizing special operators for crossover and mutation.ADP exploits the information of infeasible individuals, automatically selects, and updates the penalty factors, thus enabling a wide exploration of the search space.
In Li et al. (2019a)  The same year, Kawachi et al. (2019) introduced a DE algorithm named LSHADE equipped with an APM CHT, which is not based solely on the objective function value or the constraint violation in the search region, but also considers their balance, and applied it to the optimization of benchmark functions provided by the Congress on Evolutionary Computation (CEC) 2017.Using the right balance between constraint violation and objective function values, the proposed method outperformed a conventional APM CHT, as well as other CHTs, in most cases considered.Also, that year Datta et al. (2019) proposed the constraint handling with individual penalty (CHIP) technique, which is a combination of a common penalty-based CHT with the bi-objective EA NSGA-II.The method considers an individual penalty factor for each constraint, instead of utilizing a global factor for the overall constraint violation; nonetheless, the overall constraint violation is employed as an auxiliary objective function subject to minimization, to adaptively compute the individual penalty factors.The benefits of the proposed CHT are the capability of constraint automatic normalization and the low associated computational demand.
In Li et al. (2020) introduced the EJADE-SP algorithm for application to optimal power flow (OPF) problems, which is an enhanced adaptive DE combined with a selfadaptive penalty (SP) CHT that includes many enhancements, e.g.strategy for dynamic population reduction, crossover rate (CR) sorting mechanism, re-randomizing of CR and scale factor F parameters.The combination of this algorithm with an SP CHT achieved comparable performance to other representative algorithms such as simple DE combined with different CHTs.

Techniques based on separation of objective and constraints
In Deb (2000) developed a penalty-based CHT called feasibility rules, applicable only to population-based algorithms, e.g.genetic or other EAs, which requires no penalty parameters and is based on comparisons between pairs of individuals according to the so-called feasibility rules used in tournament selection.The major advantage of the proposed method is the elimination of the requirement of appropriate penalty factor selection by the user, which is otherwise a necessity in penalty-based CHTs to guide the search to the optimum solution.
In the same year, Runarsson and Yao (2000) introduced a novel CHT called stochastic ranking (SR) which considers the balance between objective and penalty functions via a stochastic bubble-sort algorithm and used it within an evolution strategy to solve 13 benchmark problems.Like the method proposed by Deb (2000), SR also offers the considerable advantage of eliminating the requirement of penalty factor selection by the user.
In Coello and Mezura-Montes (2002) employed within the context of structural design optimization problems the niched-pareto genetic algorithm (NPGA), which performs tournament selection based on the principle of nondomination frequently employed in multi-objective problems.For the NPGA application to single-objective COPs, the constraints are handled as additional objective functions.Using the multiobjective concept, the proposed method handles constraints efficiently without using penalty functions, while also maintaining the diversity in the population without any niching.
In Wang et al. (2009) introduced a hybrid EA that combines two mutation operators with a simplex crossover operator to generate new offsprings and employed it in 13 well-known benchmark test problems.Within the algorithm an adaptive CHT is implemented, which considers a different CHT depending on the population status based on infeasible, semi-feasible and feasible conditions.The proposed method has numerous benefits, including simplicity of implementation, enhanced robustness, and effectiveness.
In Takahama and Sakai (2010), proposed a DE algorithm with an archive and gradient-based mutation operator combined with an ε-constrained CHT, and successfully applied it to 18 test problems from CEC2010.The proposed CHT conducts ε-level comparisons between individuals with automatic control parameter adjustment, to relax the constraints and keep useful information from infeasible individuals during the search process.
In Mazhoud et al. (2013) examined a CHT called constraint violation with interval arithmetic (CVI) and employed it with a customized PSO algorithm in 24 benchmark mathematical and 3 engineering optimization problems.CVI employs the total constraint violation as an additional objective function subject to minimization, interval arithmetic for total violation normalization, as well as a simple lexicographic method for solving subsequent problems.
In Zhang et al. (2015) implemented in a new EA called backtracking search algorithm (BSA) different CHTs, such as feasibility rules and ε-constrained method with two ways of controlling the ε value and applied it to 13 benchmark functions and 4 engineering problems.The ε-constrained CHT with self-adaptive control of the ε value exhibited superior performance in respect of efficiency and convergence speed compared to the other two CHTs.
In Chehouri et al. (2016) used a GA as a numerical tool to implement a parameterfree CHT for application to NLP problems.The approach is based on the use of the violation factor and entails the population distinction into two families of feasible and infeasible solutions, where the latter contains individuals violating at least one constraint; for the population evolution at each generation, the first family is sorted with respect to the value of the fitness function, while the second according to comparison rules (feasibility-based rules).The distinctive feature of these rules is the additional consideration of the number of violated constraints, beyond the constraints violation value, in the comparison of two infeasible individuals.
In Peng et al. (2018) employed EAs combined with a CHT based on biased dynamic weights and applied it to 24 CEC2006 and 18 CEC2010 benchmark problems.Biased weights enable the selection of individuals with low objective and constraints violation values, following which the weights are dynamically adjusted to focus mostly on those individuals.In this manner, the approach prioritizes the selection of infeasible individuals from the population of parents and children, which are closer to the feasible search space.This CHT has been shown to exhibit a stable, reliable, and competitive performance when integrated in a DE algorithm.
In Fan et al. (2018) integrated an improved ε-constrained CHT called IEpsilon into DE algorithm LSHADE44 and applied it to 28 CEC2017 benchmark problems.The distinctive feature of IEpsilon compared to the classic ε-constrained handling method is the adaptive control of the ε value with respect to the proportion of feasible individuals in the current population, which balances the search between feasible and infeasible space during the evolutionary process.
In Rodrigues et al. (2018) proposed an alternative version of the balanced ranking CHT method, called Extended-BRM, and implemented it in evolutionary algorithms for application to CEC2006 and CEC2010 test functions and engineering problems.The proposed technique employs a self-adaptive mechanism to rank feasible and infeasible individuals in two separate groups, and later merges them during the search process with respect to the proportion of feasible and infeasible individuals in the current population.This CHT requires no parameter adjustment by the user, produces the most feasible solutions, and achieves superior performance compared to other CHTs, such as stochastic ranking, adaptive penalty method, etc.
In Li et al. (2019b) 2020) introduced a feasibility-rule-based CHT named FROFI, integrated into a radial basis function (RBF) surrogate-assisted DE algorithm, and applied to 2 test problems and 2 reservoir models.FROFI integrates information from the objective function through the DE algorithm operators, along with a replacement procedure, into the well-known Deb rules, to achieve a more efficient balance between objective function and constraints.
Based on the premise that among two solutions the one with better fitness function value is preferred, in Liu et al. (2020) integrated an improved Deb rule called IDeb into a Fruit fly optimization algorithm (FOA) (a search strategy based on memory, inspired by the foraging behavior of fruit flies) for application to the optimization of truss structures.The method initially employs the classic Deb rule to identify feasible individuals in the solutions memory size and subsequently evaluates the constraints violation, alas only of the individuals with better objective function values compared to the worst one stored in memory.Using the IDeb CHT a substantial reduction in the computational demand associated with structural analysis procedures can be achieved, since in such problems constraints violation evaluations are more computationally expensive compared to objective function evaluations.
In the same year, Chu et al. (2020) employed an ε-constraint CHT in a DE algorithm for topology optimization of polyline-based core sandwich structures (PCSSs).Amongst the volume, angle, and thickness constraints to be fulfilled, the predominant volume constraint is approximated by the RBF surrogate model for low computational cost, at first, and then handled through the EC method to compare different infeasible solutions; in this manner a more effective search space exploration is achieved, while potential errors arising from the prediction model RBF are reduced.Shabani et al. (2020) proposed a metaheuristic search and rescue (SAR) optimization algorithm, which simulates people behavior during search and rescue operations.A high performance was achieved by SAR when combined with an ε-constrained CHT to several benchmark test functions and 7 engineering problems.
In the work by Stanovov et al. (2020), several variations of the ε-constraint CHT are used with a DE algorithm to optimize benchmark functions and economic load dispatch problems.Among the classical ε-constraint CHT, the proposed variants of ε-constraint CHT, and the Iepsilon method, the most efficient was the combined fitness-ε vector length with individual penalty levels (EC IFL ) variant, which considers both constraints violation and objective function values, as well as each constraint individually.Rodrigues et al. (2020) proposed a PSO based algorithm combined with a weighted dynamic objective constraint handling method (WDOCHM-PSO) for design optimization of doubly fed induction generators (DFIGs).This method entails splitting the constrained problem into two unconstrained objective functions: the first considers with the distance of infeasible individuals from the feasible region (i.e.violation); and the second considers the optimization problem without constraints and is used only for feasible individuals.Mao et al. (2019) adopted the SR CHT integrated into an optimization algorithm called shuffled complex evolution (SCE), which combines the competitive complex evolution (CCE) algorithm with the simplex method to examine constrained reservoir scheduling problems.The SR CHT has been shown to be highly efficient in the identification of feasible solutions, even in the early steps of the optimization process, while its integration into SCE results in a highly robust algorithm, capable of detecting feasible regions swift.

Combination of constraint handling techniques
A popular trend in recent years is the combination of different CHTs during the optimization process.In Mallipeddi and Suganthan (2010) employed an ensemble of CHTs (feasibility rules, self-adaptive penalty, ε-constraint method, and stochastic ranking) combined with the Evolutionary Programming (EP) and DE algorithms in the context of COPs.In the proposed approach each CHT refers to a different population, alas every parent individual is compared with the produced children of all populations.A combination of different CHTs achieved superior performance compared to individual CHTs.
In Trivedi et al. (2017) adopted an ensemble of the penalty method and feasibility rules CHTs, integrated in a Unified Differential Evolution (UDE) algorithm, for the optimization of 28 CEC 2017 test problems.The proposed approach applies penalization to the first half of maximum allowed function evaluations, for efficient exploration of the search space, and feasibility rules to the other half FEs, for efficient exploitation.
In Biswas et al. (2018) examined standard IEEE 30 57-and 118-bus systems for OPF objective functions, including cost, emission, voltage stability etc., using a self-adaptive penalty CHT, a feasibility rules CHT, as well as an ensemble of these, integrated within a DE algorithm.In most cases, a combination of CHTs outperformed individual CHTs.
In the same year, Malan (2018) used a landscape-aware constraint handling approach integrated in a DE algorithm for the optimization of 18 CEC2010 test problems.The distinguishing feature of the landscape-aware approach compared to well-known ensemble of CHTs is the algorithm switch between four CHTs based on the characteristics of the landscape gathered during the search process, with only one CHT being active at a time.
In the work by Wang et al. (2019) in 2019, a DE algorithm variant called DEVNS, which produces each trial individual using two special mutation strategies and random crossover and mutation values, is combined with feasibility rules and ε-constrained CHTs in the optimization of several CEC2017 test problems.In the proposed approach, feasibility rules are applied for comparison of trial and target individuals, however all individuals with lower constraint violation than the specified ε level are considered feasible.In this manner a relaxation of constraints is achieved, and hence a more efficient search.
Also in 2019, Javed et al. ( 2019) used an ensemble of feasibility rules, self-adaptive penalty, ε-constrained and stochastic ranking CHTs, integrated into JADE and SADE algorithms, in the optimization of 24 CEC2006 benchmark problems.Both algorithms achieved very competitive feasibility rate, alas inferior success rate, compared to other considered algorithms employing a single CHT, based upon the CEC2006 evaluation criteria.
With focus on numerical optimization problems, Kaucic et al. ( 2020) introduced a PSO-based algorithm called IC-PSO, combined with a hybrid CHT which combines Deb rules with a correction mechanism.The proposed CHT utilizes a suitable operator for repairing infeasible individuals to become feasible and subsequently applies feasibility rules, thus accelerating the PSO algorithm convergence to feasible regions.

Repair algorithm-based constraint handling techniques
In the solution of COPs, a specific procedure (typically heuristic) called repair algorithm is implemented with the aim to repair an infeasible individual, i.e. convert it to a feasible one.The repaired individual can then be used solely for fitness function evaluation purposes or for replacing the initial individual into the next generation during the evolutionary process.The use of the repaired individual depends principally on the nature of the problem and the developer decision and varies between either replacing the original individual with the repaired one, or only a small percentage, or every infeasible individual in the population during the evolutionary process.
In Chootinan and Chen (2006) integrated into a simple GA a repair mechanism based on the constraint gradient information, to solve 11 optimization test problems.Constraint gradient information is obtained either via approximate finite differences or direct derivatives of the constraints, with the aim of guiding the infeasible individuals into the feasible search space.The considered CHT requires only a repair probability parameter adjustment by the user and has been shown to achieve comparable performance to other considered CHTs, such as stochastic ranking.
In the work by Salcedo-Sanz ( 2009), a comprehensive survey of the main repair heuristics used with EAs is conducted, to examine their performance in COPs.
Recently, Zein et al. (2017) employed a repair operator integrated into a GA to perform in the context of a composite structure preliminary design.The approach considers an inner and an outer loop, where the inner solves the problem using the GA and the surrogate models, while the outer reconstructs the surrogate model according to the optimal surrogate solution.The proposed method has been shown to achieve a good balance between accuracy and computational demand.
In Li et al. (2018) integrated in PSO and DE algorithms a modified repair procedure, to handle the equality constraints on Economic Dispatch problems using six different test generators systems.The repair technique evaluates the objective function derivative and uses it to adapt incrementally the output of the unit by sharing the unbalanced amount of the system constraint violation.The employed CHT has been shown to be very efficient, mainly for large-scale power systems.Gandomi and Deb (2020) in 2020 implemented a CHT called boundary update (BU) method in EAs and mathematical algorithms and tested it in several engineering and mathematical optimization problems.BU repairs the limits of selected variables to avoid constraints violation by restricting and changing them during the optimization process, generating solutions within the search space "updated" limits.Despite the BU method requiring constraint variables pre-categorizing, it achieves a high performance in the solution of COPs with complex constraints.

Boundary-based constraint handling techniques
The constraints in expression (1) can be distinguished into function constraints and bound constraints.Function constraints encompass inequality or equality functional forms of the decision variables, while bound constraints directly enforce upper and lower values (bounds) to the decision variables.While most of the research has focused on the development of CHTs for function constraint handling, a few significant research efforts have concentrated on the development of Boundary Constraint Handling Methods (BCHMs).BCHMs modify and correct the position of infeasible variable solution vectors, to reset them in the search space and convert them into feasible.
In the work by Gandomi and Yang in (2012) a comparison of an evolutionary boundary constraint handling (EBCH) scheme integrated into DE with other boundarybased methods is conducted, based upon several test problems.In the proposed method, individual components that violate a bound are replaced by a random component, which lies in between the violated bound and the corresponding component of the current best solution.The proposed scheme achieves very good performance, in most cases surpassing other classical BCHMs, such as absorbing, toroidal, reflecting, and random.
In Gandomi and Kashani (2018) proposed a probabilistic evolutionary boundary constraint handling (PEBCH) method combined with PSO for optimizing 7 benchmark test functions.PEBCH is a probabilistic version of EBCH (Gandomi and Yang 2012), which additionally considers the boundaries violation quantity.Specifically, a probability distribution function is employed, the density of which varies depending on the distance of the violated solution and boundary, taking higher values close to the violated boundary.
The results of utilizing various maps based on chaotic behavior within evolutionary algorithms (EAs) are investigated by Mozaffari et al. (2019).To handle the boundary constraints, a controlling formula which transfers a violated solution into the acceptable solution limits is used to achieve the correction of a violated solution position; in this manner global search is enhanced and early convergence is prevented.
In Biedrzycki et al. (2019) conducted a performance assessment of various combinations between seven DE algorithms and seventeen BCHMs, based upon CEC2017 benchmark functions, concluding that the efficiency of DE algorithms varies substantially depending on the employed BCHM.
In the same year, Juárez-Castillo et al. ( 2019) implemented an adaptive BCH scheme in PSO and DE algorithms for the optimization of sixty single-objective COPs.The adaptive approach employs a specific BCHM if the population consists entirely of infeasible solutions, otherwise a BCHM from a predefined set is randomly selected, based on a dynamically-updated probability associated with it.This adaptive scheme achieved better performance compared to single BCHMs in COPs with moderate or large number of variables or constraints, as well as in COPs where the identification of feasible individuals is rendered cumbersome.
In Arouri and Sayyafzadeh (2020) developed a gradient-based algorithm that combines features of simultaneous perturbation stochastic approximation (SPSA) and an adaptive moment estimation framework, for application to production optimization problems.The algorithm was tested using the projection and logarithmic transformation BCHMs, with the former achieving superior performance compared to the latter in both cases examined.
In the work by Biedrzycki et al. in the same year (2020), the performance of the covariance matrix adaptation evolution strategy (CMA-ES) algorithm is assessed in the optimization of unimodal and multimodal CEC2017 and Black-Box Optimization Benchmarking (BBOB) problems, considering 22 BCHMs categorized as repair, feasibility preserving, specialized and penalty functions.Amongst all considered BCHMs, Darwinian reflection and resampling achieved the highest performance.

Other constraint handling techniques
In Atamna et al. (2020) implemented an augmented Lagrangian approach as CHT in Evolution Strategies for optimizing sphere and ellipsoid functions with linear constraints, where the initial constrained objective function is converted into an unconstrained augmented Lagrangian function, while its parameters are adopted during optimization.
In the same year, Qian et al. ( 2020) employed a GA combined with a surrogate model for solving 9 numerical and 2 engineering optimization problems, using the Kriging surrogate prediction model instead of the actual constraints, which is continuously updated during the optimization procedure to avoid infeasible optimization solutions.This method leads to optimal feasible solutions and is associated with a low computational demand compared to other existing techniques.Rosso et al. (2021) presented a new nonpenalty-based constraint handling approach for PSO type of algorithms adopting support vector machine (SVM), that is a supervised classification machine learning approach.Due to its generality, constraint handling with SVM appears more adaptive both to nonlinear and discontinuous boundary, while aiming to preserve the feasibility of the population, a simple bisection algorithm was also implemented.
While recently, Rosso et al. (2022) a variant of the well-known swarm-based algorithm, the Particle Swarm Optimization (PSO), is developed to solve constrained problems with a different approach to the classical penalty function technique.Stateof-art improvements and suggestions are also adopted in the current implementation (inertia weight, neighbourhood).Furthermore, a new local search operator has been implemented to help localize the feasible region in challenging optimization problems.This operator is based on hybridization with another milestone meta-heuristic algorithm, the Evolutionary Strategy (ES).The self-adaptive variant has been adopted because of its advantage of not requiring any other arbitrary parameter to be tuned.

Constraint handling technique formulation
In this section, four of the most commonly employed CHTs identified in the literature are described, specifically: (i) penalty methods; (ii) feasibility rules; (iii) ε-constrained method; and (iv) stochastic ranking, while a brief description of methods based on ensemble of CHTs (i)-(iv) is also provided.Furthermore, eight novel variants of the feasibility rules and ε-constrained CHTs are introduced.The performance of these novel variants is evaluated later in Sect.5, based upon 20 mathematical and engineering COPs, and compared to four existing variants of CHTs (i)-(iv), specifically: (i) adaptive penalty by Kawachi et al. (2019) (ii) feasibility rules by Deb (2000); (iii) improved ε-constrained by Fan et al. (2018); and (iv) stochastic ranking by Runarsson and Yao (2000), which are also described herein.

Penalty methods
In penalty methods the constrained problem is converted into an unconstrained one with the addition of a penalty term to the objective function.The penalty term depends on the level of violation of the constraint functions, as well as penalty factors whose value can vary throughout the optimization process.The penalized objective function can be expressed in the following general form (Coello and Montes 2002): where F(x) is the penalized objective function (also called fitness function), g j (x) and h k (x) are the inequality and equality constraints, respectively, and r j and c k are positive weight coefficients called penalty factors.Converting the equality into inequality constraints, as previously discussed in Sect. 1, Eq. ( 2) can be rewritten as: (3) Kawachi et al. (2019) proposed an adaptive procedure for calculating the penalty factor during the evolution process, which effectively deals with issues of over-or under-penalization that potentially led to divergence of the search process from the optimum.More specifically, a fitness function is formulated: where υ(x) is the mean constraint violation and P F is the penalty factor, and subsequently three steps are implemented.First the penalty factor candidates (PFCs) are calculated by comparing two individuals as follows: where subscripts k and l denote two distinct individuals; P FC is calculated for all possible combinations in a given population.In the second step, the penalty factor (P F) is defined, as follows: if the percentage of negative P FCs exceeds 50%, the P F remains the same with the one of the previous generations; otherwise P F is calculated as the average value of the positive P FCs.In the third step P F is updated.
If the proportion of feasible individuals r f in the population exceeds the p f eas , then the penalty factor of the next generation is obtained as follows: where p rate ∈ [0, 1] and p f eas are user defined parameters.

Feasibility rules
In penalty methods, different values of the penalty factor need to be tested to establish the best option, since inappropriate values can cause search process divergence from the area of the optimum solution.To overcome this shortcoming Deb (2000) proposed a new technique based upon a tournament selection operator for comparing two solutions, in accordance with the following criteria: (i) any feasible solution is preferred to any infeasible solution; (ii) among two feasible solutions, the one associated with better objective function value is preferred; and (iii) among two infeasible solutions, the one associated with smaller constraint violation value is preferred.

Original feasibility rules technique
According to the original feasibility rules technique, the fitness function is defined as follows: where f max is the objective function value for the worst feasible solution of the current population.If no feasible solution exists in a population, f max is set to zero.According to this technique, all constraint violations are summed, and this sum is compared as a single value.Therefore, in case of infeasible solutions, these are compared based solely on the level of constraint violation.

New variants of the feasibility rules technique
In this section four new variants of the feasibility rules technique are presented.Specifically, the comparison rules described above for the original technique have been modified to bias the search direction towards better feasible solutions, instead of simply feasible ones, to achieve a global exploration of the design space.The formulation of the proposed variants relies on the premise that a small constraint violation should not be the sole criterion of selecting one individual over another, as is the case in existing techniques.Specifically, when comparing two infeasible individuals, the objective function value and the number of violated constraint functions should also be considered, alongside the maximum constraint violation value, and the product of these three entities shall determine the final winner: if an individual is associated with a lower level of constraint violation than another, but the objective function value of the first is much larger compared to the second, maintaining the first one in the next generation might direct the search process to a local optimum.Of course, if the second individual corresponds to a very large level of constraint violation, it will be penalized with a much larger penalization value.To avoid very small objective function values that could misdirect the search process, the objective function value of the best feasible individual found so far is used as the lower bound.Analogously, when comparing one feasible and one infeasible individual, the second one could be chosen based on the product of its objective function value and the level of the constraint violation.If the objective function value of the feasible solution is very bad, then selecting it might direct the search process to a local optimum, and hence the infeasible individual should be preferred.Finally, since it is not sure which one is better between a solution with many violated constraint functions with small level of violation compared to a solution with just one violated constraint function of very large level of violation, the total number of violated constraint functions is also considered together with the maximum level of constraint violation.Taking the above into consideration, the formulation of the new variants is described hereafter.Please note that the red colour value has been changed to italic in Tables 8, 15, 18, 21, 22. Please check and confirm.Confirm To calculate the fitness function of an infeasible individual, the p violation factor is introduced.This is the normalized maximum constraint violation of the individual, which, depending on the variant, is multiplied with a term related to the relative number of violated constraint functions for a solution over the total number of constraint functions.The p violation factor for each of the four variants proposed is defined as follows: where g j (x) is the normalised value of the j th constraint (an example of a normalised constraint is shown below in Eq. 9b), n constviol is the number of violated constraint functions and n const is the total number of constraint functions.The fitness function of an individual is calculated as the product of p violation with the maximum among the objective function value of the individual and the best objective function value of the feasible individuals found so far, as stated in the following expression: where an example of a constraint function is shown in Eq. (9b), where σ (x) is the principal stress developed into a specific structural element and σ all is the corresponding allowable value.On the basis of expression (9b), three criteria are used to compare two solutions: (i) between two feasible solutions, the one associated with better objective function value is preferred, (ii) between two infeasible solutions, the one associated with lower level of constraint violation and lower number of violated constraint functions is preferred; and (iii) between one feasible and one infeasible solution, the selection of the feasible solution depends on its objective function value, the level of constraint violation and the number of violated constraint functions of the infeasible one.In case of a maximization problem the criteria are adapted accordingly.

ε-constrained method
A modified version of the feasibility rules technique was originally proposed by Takahama and Sakai (2005), with the aim of providing comparison rules on constrained optimization problems during the evolutionary process.According to this technique, the overall constraint violation function ϕ(x) is defined as either the maximum or the sum of all constraints, and is used in conjunction with the ε level comparison {< ε }; for any ε ≥ 0, an individual ( f 1 , ϕ 1 ) is considered better than another individual ( f 2 , ϕ 2 ) according to the following comparison rules: Therefore: (a) if both overall constraint violations ϕ(x) of two individuals are either less than or equal to ε, the one with better objective function value f is preferred; (b) If the overall constraint violations ϕ(x) of two individuals are equal, the one with better objective function value f is preferred; and (c) if either or both constraint violations ϕ(x) of two individuals are larger than ε, the one with lower constraint violation ϕ(x) is preferred.Building upon their previous work, Takahama and Sakai (2010) suggested that the ε level can be usually controlled according to the following expression: where the initial ε level is equal to ϕ θ , that is the constraint violation of the top θ th individual in the initial population, t is the number of iterations, T C is the control generation and c p is a user defined parameter with value greater or equal to 3.

Improved ε-constrained method
In Fan et al. (2018) proposed the following improved ε setting approach, which applies the ε level comparison {< ε } to establish comparison rules during the evolutionary process, and further balance the evolutionary search of the population between feasible and infeasible regions: where r k is the proportion of feasible solutions (PFS) in the generation; ϕ θ is the θ th largest overall constraint violation of all individuals in the initial population (sorted in an array by their overall constraint violation), where θ = γ • N pop , N pop being the population size and γ ∈ [0.2, 0.8]; ϕ max is the largest overall constraint violation of all individuals; T C is the termination function evaluations (FEs), where T C ∈ [0.1max F Es, 0.8max F Es]; and c p ∈ [2, 10], τ ∈ [0, 1] and α ∈ [0, 1] are user-defined parameters.

New variants of the ε-Constrained method
Analogously to the new feasibility rules variants introduced previously in Sect.3.2.2,new variants of the ε-constrained technique are also proposed herein.The ε level comparison {< ε } in the new variants for two individuals 1, 2 is defined via the following four rules: where f b f stands for the objective function value of the best feasible individual found so far, and p 1 , p 2 denote the value of the p violation factor calculated for individuals 1 and 2, based upon Eqs.(8a) to (8d).Therefore: (a) if both maximum constraint violations of two individuals are either less than or equal to ε level, the one associated with better objective function value f is preferred; (b) if the maximum constraint violations of two individuals are equal, the one associated with better objective function value f is preferred; (c) if both maximum constraint violations of two individuals are larger than the ε level, the one associated with better product of constraint violation and objective function value is preferred; and (d) if only one maximum constraint violation of two individuals is larger than the ε level, the other individual is preferred only if its objective function value is better than the product of the constraint violation with the objective function value of the second individual.In line with the new variants of the feasibility rules technique, in addition to the maximum constraint violation value the total number of violated constraint functions is also considered.To control the ε level, the formulation proposed by Fan et al. ( 2018) is adopted, as given in Eq. ( 12).

Stochastic ranking method
A frequently employed CHT that circumvents the shortcomings of penalty techniques is the stochastic ranking method, originally proposed by Runarsson and Yao (2000).Stochastic ranking is based upon comparing all adjacent individuals in a population according to penalty and objective function dominance, to rank them from best to worst.
A probability P f ∈ [0, 1] associated with the sole use of the objective function for the adjacent individual comparisons must be defined by the designer to achieve enhanced performance.A value P f = 1/2 means that the probability of a comparison based on the objective function is the same as the probability based on the penalty function.Ranking is achieved via a bubble-sort-type algorithm, where λ individuals are ranked by comparing adjacent individuals in N sweeps (the least λ sweeps) in every generation.If both adjacent individuals are feasible or a randomly generated number u is less than the selected probability P f , the individuals are compared according to their objective function values; otherwise, if one or both are infeasible, the comparison is based on the constraint violation value.The algorithm is terminated if no change takes place in the ranking order within a complete sweep.

Ensemble of constraint handling techniques
It is not always easy to determine which CHT is better for determining the global optimum in a given COP.While a specific CHT might be more suitable than others for a particular stage of the evolution process, it might be inappropriate for another one, and hence different CHTs might be more effective during different stages of the search process.This depends on several factors, such as the multimodality of the problem, the proportion of feasible search region in the whole search region, the adopted algorithm etc. Motivated by these observations, Mallipeddi and Suganthan (2010) proposed a combination of four CHTs (feasibility rules, self-adaptive penalty, ε-constraint method, and stochastic ranking), a scheme called ECHT, to examine constrained optimization problems, where each technique handles a different population.Parents and offsprings of the population compete not only among each other, but also with the respective entities of all the other populations.This means that if an offspring generated from a particular population is not selected for the next generation based on this population's specific CHT, it could be selected by other populations that use a different CHT.

Optimization computing platform and pity beetle algorithm
The solution of real-world, large-scale optimization problems within the field of structural engineering can only be achieved with a synergy of the following aspects (Melchiorre et al. 2021;Lagaros et al. 2022;Cucuzza et al. 2022): (i) efficient numerical modelling of the physical problem; (ii) reliable numerical optimization algorithms for enhanced structural design; (iii) rationalized modelling of the system uncertainties (in the case of probabilistic design problems); and (iv) exploitation of modern High-Performance Computing (HPC) facilities.The original OCP (Lagaros 2014) and its evolution HP-OCP offer all aforementioned capabilities, including efficient numerical modelling tools and several MOAs with a great variety of embedded features and CHTs, while its implementation enables its straightforward upgrade to incorporate new tools and features.
In Kallioras et al. (2018) introduced the novel nature-inspired pity beetle MOA, a single-objective optimization algorithm for unconstrained problems which has been shown to handle NP-hard optimization problems efficiently.PBA is inspired by the behavior of the six-toothed bark beetle (pityogenes chalcographus beetle) when searching for food on the bark of trees and simulates the beetle behavior in three stages: (i) population initialization; (ii) population regeneration; and (iii) location update.In the initialization stage (i), a population consisting of males and females is randomly allocated in the search space via a random sampling technique (RST), which ensures adequate diversification to avoid premature convergences.The male particles act as pioneers in search for the most suitable host, producing pheromone that attracts the other males and females.In the population regeneration stage (ii), each single particle seeks an improved position in the search space to create its own population.This is achieved via five types of hypervolume selection patterns, namely global-scale, large-scale, mid-scale, neighboring search, and memory consideration search volume patterns, where in the latter the best positions identified are stored in memory.In the location update stage (iii), the position of each mating male and female is updated and the previous positions are eliminated, except for those stored in memory (memory consideration search volume).Collectively, the hypervolume patterns are selected in accordance with the following expression: In expression ( 14): g refers to the current generation (pursuit step); i signifies the component of the design vector; j denotes the member of the population; k designates the brood (colony); n denotes the number of unknowns; N P O P is the population size; N broods is the maximum number of broods (termination criterion); x (g) birth,i is the i th element of the birth position/solution vector; x (g) j,k is the new position vector found in population k; superscripts l, u denote the lower and upper bound; while the algorithmic parameters of PBA are: the population of pioneer particles N pop ; the neighboring factor f nb ; the fine tuning factor f tn , which is used to define the area size where a very narrow local search is performed, in case a better position is identified; the mid-scale factor f ms ; the large-scale factor f ls ; the probability of choosing large-scale search or memory consideration pr ; and the function evaluations multiplication factor f F E , which is used to define the maximum number of unsuccessful function evaluations F E un with reference to the total function evaluations F E un = F E total • f F E ; and MEM denotes the memory of size N pop , where the best solutions found by the algorithm up until the current step are stored in memory.A detailed description of PBA can be found in Kallioras et al. (2018).

Numerical examples
In this section, the performance of the 8 novel feasibility rules and ε-constrained CHT variants proposed in Sect. 3 is evaluated and compared against 4 well-known CHT variants from the literature: (i) adaptive penalty by Kawachi et al. (2019) (ii) feasibility rules by Deb (2000); (iii) improved ε-constrained by Fan et al. (2018); and (iv) stochastic ranking by Runarsson and Yao (2000).Consideration is given to 20 COPs previously investigated by other authors, encompassing: (i) 14 mathematical COPs, including a 2D problem investigated by Deb (2000) and 13 problems investigated by Runarsson and Yao (2000) (labelled G1-G13); and (ii) 6 engineering COPs, including a welded beam, pressure vessel and tension-compression string design investigated by Tsipianitis and Tsompanakis (2020), as well as 2D and 3D truss design examples investigated by Farshi and Alinia-ziazi (2010).
For each numerical example and each employed CHT, 20 independent optimization computations are performed (20 COPs × 12 CHTs × 20 computations = 4800 total computations), to establish an adequate sample for a statistical analysis of the obtained results.In all computations the PBA (Kallioras et al. 2018) implemented in the MATLAB variant of HP-OCP is employed (Sect.4), using identical parameters, termination criteria and bound-constraint treatment (out of bounds solution vectors corrected via projection to the closest bound), which eliminates the influence of the PBA specifications, thus enabling a common basis for CHT comparison to be established.The adopted PBA parameters are: N pop = 30, f nb = 0.08, f tn = 0.015, f ms = 0.90, f ls = 100, f F E = 0.25 and pr = 0.20 (Kallioras et al. 2018), as well as F E un = 20, 000.The specifications of the employed CHT variants are given in Table 2.All computations were conducted using MATLAB ® version R2020b, on an x64 PC with 14-Core Intel ® Xeon ® processor with 128 GB RAM memory.

Mathematical constrained optimization problems
The first mathematical COP is a 2D minimization problem originally introduced by Deb (2000), which is formulated as follows: subject to : The reported objective function value in Deb (2000) is equal to 13.59658.The results obtained using the PBA and the various CHTs are given in Table 3.The best performance was achieved by AdPenalty, EconstOR, FeasRulesNEW4 and Econst-NEW2 CHTs with an optimized objective function value lower than 13.60, while EconstNEW2 achieved the lowest CoV value of 0.142% Subsequently, consideration is given to 13 mathematical COPs labelled G1-G13, originally introduced by Runarsson and Yao (2000).The problem specifications are given in Table 4, where LE, LI, NE and NI respectively denote the number of Linear Equality, Linear Inequality, Nonlinear Equality and Nonlinear Inequality constraint functions.The results obtained using the PBA and the various CHTs are given in Tables 5,6,7,8,9,10,11,12,13,14,15,16 and 17.G1 is a 13-dimensional minimization problem with 9 constraint functions and reference objective function value equal to − 15.0 (Runarsson and Yao 2000).The respective results are given in Table 5.All CHTs produced optimized objective function values below − 14.0 and the highest performance was achieved by FeasRulesNEW4 with − 14.99.The lowest CoV value of 1.36% was achieved by EconstNEW4.G2 is a 20-dimensional maximization problem with 2 constraint functions and reference objective function value equal to 803,619 (Runarsson and Yao 2000).The respective results are given in Table 6.Almost half of the CHTs produced optimized objective function values greater than 0.75 and the highest performance was achieved by EconstNEW2 with 0.7992425.The lowest CoV value of 1.74% was achieved by StochRanking.
G3 is a 20-dimensional maximization problem with 1 constraint function and reference objective function value equal to 1.0 (Runarsson and Yao 2000).The respective results are given in Table 7. Almost all CHTs produced optimized objective function values greater than 0.95 and the highest performance was achieved by EconstNEW2 with 0.995233.The lowest CoV value of 0.66% was achieved by EconstNEW4.
G4 is a 5-dimensional minimization problem with 6 constraint functions and reference objective function value equal to − 30,665.539(Runarsson and Yao 2000).The respective results are given in Table 8. Almost all CHTs produced optimized objective function values lower than − 30,000 and the highest performance was achieved by AdPenalty with − 30,541.540.Nearly all techniques achieved CoV values less than 1.0%, with the lowest corresponding to EconstNEW1, StochRanking failed to converge.
G5 is a 4-dimensional minimization problem with 5 constraint functions and reference objective function value equal to 5126.4981 (Runarsson and Yao 2000).The respective results are given in Table 9. Almost all CHTs produced optimized objective function values lower than 5300 and the highest performance was achieved by Econst-NEW3 with 5127.392.The lowest CoV value of 0.98% was achieved by EconstNEW1.
G6 is a 2-dimensional minimization problem with 2 constraint functions and reference objective function value equal to − 6961.81388(Runarsson and Yao 2000).The respective results are given in Table 10.Almost half CHTs produced optimized objective function values lower than − 6700 and the highest performance was achieved by FeasRulesNEW1 with − 6912.992.The lowest CoV values of 3.85% and 4.01% were respectively achieved by FeasRulesOR and EconstNEW2.
G7 is a 10-dimensional minimization problem with 8 constraint functions and reference objective function value equal to 24.3062091 (Runarsson and Yao 2000).The respective results are given in Table 11.Almost all CHTs produced optimized objective function values lower than 25 and the highest performance was achieved by FeasRule-sOR with 24.335.Nearly all techniques achieved CoV values less than 4.0%; the lowest CoV values of 2.03% and 2.16% were respectively achieved by EconstNEW4 and FeasRulesNEW2.
G8 is a 2-dimensional maximization problem with 2 constraint functions and reference objective function value equal to 0.095825 (Runarsson and Yao 2000).The respective results are given in Table 12.Almost all CHTs produced optimized objective function values greater than 0.090 and the highest performance was achieved by EconstOR with 24.335.The lowest CoV values of 0.13% was also achieved by EconstOR.
G9 is a 7-dimensional minimization problem with 4 constraint functions and reference objective function value equal to 680.6300573 (Runarsson and Yao 2000).The respective results are given in Table 13.Almost all CHTs produced optimized objective function values greater than 750 and the highest performance was achieved by FeasRulesNEW3 with 738.137.The lowest CoV value of 1.50% was achieved by EconstOR.
G10 is a 8-dimensional minimization problem with 6 constraint functions and reference objective function value equal to 7049.3307 (Runarsson and Yao 2000).The respective results are given in Table 14.Almost all CHTs produced optimized objective function values lower than 7400 and the highest performance was achieved by EconstNEW1 with 7058.463.The lowest CoV value of 1.07% was achieved by FeasRulesNEW2.
G11 is a 2-dimensional minimization problem with 1 constraint functions and reference objective function value equal to 0.75 (Runarsson and Yao 2000).The respective results are given in Table 15.Almost all CHTs produced optimized objective function values lower than 0.78 and the highest performance was achieved by AdPenalty, EconstOR, EconstNEW1, FeasRulesNEW2, EconstNEW2 and EconstNEW4 with values very close to 0.75.The lowest CoV value of 3.11% was achieved by AdPenalty.StochRanking failed to converge.
G12 is a 3-dimensional minimization problem with 729 constraint functions and reference objective function value equal to 1 (Runarsson and Yao 2000).The respective results are given in Table 16.All CHTs produced the reference objective function value, with most of them being associated with an extremely low CoV below 0.02%.
G13 is a 5-dimensional minimization problem with 3 constraint functions and reference objective function value equal to 0.0539498 (Runarsson and Yao 2000).All CHTs produced optimized objective function values lower than 0.056 and the highest performance was achieved by FeasRulesOR, EconstNEW1 and EconstNEW4 with values very close to 0.054.The lowest CoV value of 1.78% was achieved by EconstNEW2.

Welded beam design
This COP concerns the optimal design of a welded beam via minimization of the fabrication cost, which depends on four design variables: weld thickness x 1 , weld length x 2 , beam height x 3 , and beam width x 4 .The problem is formulated as follows: Minimi ze f (x) = 1.10471x 2 1 x 2 + 0.04811x 3 x 4 (14.0 + x 2 ), subject : where . Identical parameters with (Tsipianitis and Tsompanakis 2020) are adopted, with P = 6000 lb, L = 14 ∈, E = 30 × 10 6 psi, G = 12 × 10 6 psi, τ max = 13600 psi, σ max = 30000 psi, δ max = 0.25 ∈.The reference objective function value is equal to 1.72485084 lb.The results obtained using the PBA and the various CHTs are given in Table 18, where evidently the best solution was achieved by EconstNEW3, with an objective function value lower of 1.729478, and the lowest CoV by AdPenalty, Econ-stOR and EconstNEW2, with values below 2%.It should be noted that StochRanking failed to converge to a feasible design solution.

Tension-compression string design problem
This numerical example concerns the weight minimization of a tension-compression string, which depends on three design variables: string diameter x 1 , mean coil diameter x 2 , and coil number x 3 .The COP is formulated as follows: Identical parameters with (Tsipianitis and Tsompanakis 2020) are adopted.The reference objective function value is equal to 0.012665 lb.The results obtained using the PBA and the various CHTs are given in Table 20.Almost all CHTs lead to weight objective function values less than 0.013 lb, with the lowest value of 0.0126911 lb achieved by FeasRulesNEW1, while CoV values less than 2.0% were achieved by FeasRulesNEW4 and EconstNEW4.

2D and 3D truss design
Three sizing COPs previously investigated in Farshi and Alinia-ziazi (2010) are examined herein, concerned with the optimal structural design in respect of weight of the 2D 10-bar, 3D 25-bar and 3D 72-bar steel truss structures shown in Figs. 1, 2 and 3.The design variables are the cross-sectional areas of individual structural elements or element groups.
For the 2D 10-bar truss problem shown in Fig. 1 an independent design variable is employed for each bar, continuously defined in the interval [0.1, 33.5] in 2 .The COP is formulated with 28 constraint functions, which impose: (i) tensile and compressive stress limits of 25 ksi for each bar; and (ii) horizontal and vertical displacement limits of 2.0 at each unrestrained node.More details on the problem formulation can be found in Farshi and Alinia-ziazi (2010).
The reference structural weight objective function value found in the literature is equal to 5057.88 lb (Farshi and Alinia-ziazi 2010).The results obtained using the PBA and the various CHTs are given in Table 21.Evidently, most CHTs achieved objective function values less than 5500 lb, with AdPenalty achieving the best solution of 5293.971lb.FeasRulesOR, FeasRulesNEW1, EconstNEW1, FeasRulesNEW3 and FeasRulesNEW4 achieved the lowest CoV, with values below 3.0%.StochRanking failed to converge to a feasible design.
For the 3D 25-bar truss problem shown in Fig. 2a-d, the structural elements are allocated to 8 groups, each corresponding to a design variable continuously defined in in the interval [0.01, 3.4] in 2 .The CPO is formulated with 110 constraint functions, which impose: (i) a tensile stress limit of 35 ksi and a compressive stress limit in accordance with the AISC code; and (ii) a limit of 0.35 on the absolute value of the displacement components of each unrestrained node.More details on the problem formulation can be found in Farshi and Alinia-ziazi (2010).The reference structural weight objective function value found in the literature is equal to 545.175 lb (Farshi and Alinia-ziazi 2010).The results obtained using the PBA and the various CHTs are given in Table 22.Evidently, most CHTs achieved objective function values less than 560 lb.FeasRulesNEW2, FeasRulesNEW4 and EconstNEW4 achieved the lowest CoV, with values below 1.0%.StochRanking failed to converge to a feasible design.
For the 72-bar truss problem shown in Fig. 3a-d, the structural elements are allocated to 16 groups, each corresponding to a design variable continuously defined in the interval [0.1, 3.0] in 2 .The CPO is formulated with 264 constraint functions which impose: (i) tensile and compressive stress limits of 25 ksi for each bar; and (ii) a limit of 0.35 on the absolute value of the displacement components of each unrestrained node.More details on the problem formulation can be found in Farshi and Alinia-ziazi (2010).The results obtained using the PBA and the various CHTs are given in Table 22.
The reference structural weight objective function value found in the literature is equal to 379.66 lb (Farshi and Alinia-ziazi 2010).The results obtained using the PBA and the various CHTs are given in Table 23.FeasRulesNEW2 and EconstNEW2 achieved objective function values less than 400 lb.AdPenalty achieved the lowest CoV of 0.54%.

Discussion
Evidently, the novel CHT variants overall outperform the original ones in most COPs, in respect of both objective function values and CoV.This is attributed to their novel feature of performing individual comparisons based not only upon the objective function value, but also the level of constraint violation and the number of violated constraints.At this point it must be mentioned that the derived conclusions concern the 20 numerical examples examined herein and the specific PBA and CHT variant implementation and can thus not be generalized.

Conclusions
This paper presents a comprehensive state-of-the-art review on almost 60 existing CHT variants that can be implemented in MOAs to solve single-objective COPs.Furthermore, eight new variants of the well-known feasibility rules and ε-constrained CHTs are proposed.The distinctive feature of the proposed variants compared to the original ones is that necessary comparisons between individuals depend not only on the objective function value, but also on the level of constraint violation and the number of violated constraint functions.The proposed variants were successfully implemented in the metaheuristic PBA (Kallioras et al. 2018), within HP-OCP, and their performance was evaluated and compared to that of the well-known adaptive penalty (Kawachi et al. 2019), original feasibility rules (Deb 2000), improved ε-constrained (Fan et al. 2018) and original stochastic ranking (Runarsson and Yao 2000) variants, based upon 20 single-objective benchmark mathematical and engineering COPs.The results clearly demonstrate the performance superiority of the proposed novel variants in comparison with the existing ones in most COPs, in respect of both objective function values and CoV values.This highlights the effectiveness and versatility of the proposed approach that considers the level and number of constraint violations for selecting individuals.Concluding, it is important to emphasize that there is no unique constraint handling technique that consistently performs best in any type of problem.
employed a Kriging model in a dynamic surrogate-based optimization (DSBO) of standard test functions and an engineering problem, utilizing two criteria for sample selection: maximization of the expected improvement (EI) function; and minimization of surrogate model prediction, and achieved constraint handling in three different ways: (a) constraining EI function, (b) penalizing the prediction of the surrogate model and (c) penalizing the objective function.DSBO achieved best performance in two tests and a beam optimization problem with the simultaneous employment of methods (a) and (b).

Fig. 1
Fig. 1 Geometric, boundary condition, material and loading specifications for 2D 10-bar truss constrained optimization problem

Fig. 2
Fig. 2 Geometric, boundary condition, material and loading specifications for 3D 25-bar truss constrained optimization problem: a perspective view; b elevation on X-Z plane; c elevation on Y-Z plane; and d plan view on X-Y plane

Table 1
Overview and classification of constraint handling techniques (CHTs) for metaheuristic optimization algorithms

Table 2
Specifications of constrained handling techniques employed in numerical examples

Table 3
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem by Deb

Table 4
Specifications of mathematical constrained optimization problems G1-G13

Table 5
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G1

Table 6
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G2

Table 7
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G3

Table 8
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G4 (solutions denoted in italic correspond to infeasible designs)

Table 9
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G5

Table 10
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G6

Table 11
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G7

Table 12
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G8

Table 13
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G9

Table 14
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G10

Table 15
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G11 (solutions denoted in italic correspond to infeasible designs)

Table 16
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G12

Table 17
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in optimization problem G13

Table 18
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in the welded beam design optimization problem (solutions denoted in italic correspond to infeasible designs)

Table 19
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in the pressure vessel design optimization problem

Table 20
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in the tension-compression string design optimization problem

Table 21
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in the 2D 10-bar truss design optimization problem (solutions denoted in italic correspond to infeasible designs)

Table 22
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in the 3D 25-bar truss design optimization problem (solutions denoted in italic correspond to infeasible designs)

Table 23
Best, median, worst, and mean objective function values and coefficient of variation obtained for the various constraint handling techniques in the 3D 72-bar truss design optimization problem