Efficiently handling constraints in mixed-integer nonlinear programming problems using gradient-based repair differential evolution

Differential evolution

Storn & Price (1997) proposed a population-based optimization algorithm, DE, that has gained popularity for its high success rate and straightforward implementation. DE consists of four fundamental processes: initialization, mutation, crossover, and selection, which are repeated in each cycle.

The initial population P is created through random sampling from the search space. Each individual in the population, denoted as ${\mathbf{x}}_i=(x_{i,1},...,x_{i,j})$, represents a solution to the problem with $x_i$ being the $i$th individual, $i=1,...,\mathrm{N}\mathrm{P}$. $j$ is the number of decision variables and NP is the population size. A mutant vector ${\mathbf{v}}_i=(v_{i,1},...,v_{i,j})$ is generated for every individual ${\mathbf{x}}_i$ as the base of the evolution. Several mutation operators have been proposed (Abbasa, Ahmadb & Jabeenc, 2017). In this study, the following versions were employed:

• DE/rand/1:

  (2) $${\mathbf{v}}_i={\mathbf{x}}_{r1}+F({\mathbf{x}}_{r2}-{\mathbf{x}}_{r3})$$

• DE/current-to-rand/1:

  (3) $${\mathbf{v}}_i={\mathbf{x}}_i+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}({\mathbf{x}}_{r1}-{\mathbf{x}}_i)+F({\mathbf{x}}_{r2}-{\mathbf{x}}_{r3})$$

• DE/best/1:

  (4) $${\mathbf{v}}_i={\mathbf{x}}_{best}+F({\mathbf{x}}_{r1}-{\mathbf{x}}_{r2})$$

• DE/rand-to-best/1:

  (5) $${\mathbf{v}}_i={\mathbf{x}}_{r1}+\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}({\mathbf{x}}_{best}-{\mathbf{x}}_{r2})+F({\mathbf{x}}_{r3}-{\mathbf{x}}_{r4})$$

DE executes a crossover operation between the target vector ${\mathbf{x}}_i$ and its mutant vector ${\mathbf{v}}_i$ to generate the trial vector ${\mathbf{u}}_i=(u_{i,1},...,u_{i,j})$. The binomial crossover is defined in Eq. (6), where the control crossover parameter CR is a value in $[0,1]$. A uniformly distributed random number rand_j is generated for each decision variable $j$ in the interval $[0,1]$, and $j_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}}$ is a random integer in the range $[1,\mathrm{N}\mathrm{P}]$.

(6) $$u_{i,j}=\{\begin{array}{cc}{v}_{i,j},\hfill & \mathrm{i}\mathrm{f}\mathrm{r}\mathrm{a}\mathrm{n}{\mathrm{d}}_j\le \mathrm{C}\mathrm{R}\mathrm{o}\mathrm{r}j=j_{\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}.}\\ x_{i,j},& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\end{array}$$

During the selection step, the best individual is chosen between x_i and u_i to be carried forward to the next generation. A set of feasibility rules (Deb, 2000) is employed to compare them as follows:

• 1. Between two infeasible solutions, the one with lower constraint violation is preferred.

• 2. If one solution is infeasible and the other one is feasible, the feasible solution is preferred.

• 3. Between two feasible solutions, the one with better objective function value is preferred.

Gradient-based repair method

The technique of gradient-based repair was first introduced for genetic algorithms (Chootinan & Chen, 2006). Takahama & Sakai (2009) incorporated this approach in a modified proposal version of the constraint handler $ϵ$-constrained. This method uses the gradient information of the constraint set to guide the infeasible solutions toward feasible regions. The gradient can be obtained directly from the constraints or, if the analytical derivation is challenging, it can be estimated by numerical techniques such as the finite differences (Chapra & Canale, 2010).

$\mathbf{V}(\mathbf{x})$ is defined by Eq. (7) as a vector that contains the violation degree for each inequality ( $g$) and equality ( $h$) constraint in a given problem, for a particular solution vector x. Parameters $n$ and $m$ are the number of inequality and equality constraints, respectively, and $ϵ$ specifies the tolerance for equality constraints. Preserving the sign of the equality violation is a crucial aspect of the method and is represented by the function $\mathrm{s}\mathrm{g}\mathrm{n}(h)$.

(7) $$\mathbf{V}\mathbf{(}\mathbf{x}\mathbf{)}=\left[\begin{array}{c}max(g_1(\mathbf{x}),0)\\ ⋮\\ max(g_n(\mathbf{x}),0)\\ \mathrm{s}\mathrm{g}\mathrm{n}(h_1)\cdot max(|h_1(\mathbf{x})|-ϵ,0)\\ ⋮\\ \mathrm{s}\mathrm{g}\mathrm{n}(h_m)\cdot max(|h_m(\mathbf{x})|-ϵ,0)\end{array}\right].$$

The gradient matrix of these constraints with respect to the N components of x, denoted as $\mathrm{\nabla }\mathbf{V}(\mathbf{x})$, is defined by Eq. (8):

(8) $$\mathrm{\nabla }\mathbf{V}(\mathbf{x})=\left[\begin{array}{cccc}\frac{\mathrm{\partial }{g}_1(\mathbf{x})}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{g}_1(\mathbf{x})}{\mathrm{\partial }{x}_2}& \dots & \frac{\mathrm{\partial }{g}_1(\mathbf{x})}{\mathrm{\partial }{x}_N}\\ ⋮& ⋮& & ⋮\\ \frac{\mathrm{\partial }{g}_n(\mathbf{x})}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{g}_n(\mathbf{x})}{\mathrm{\partial }{x}_2}& \dots & \frac{\mathrm{\partial }{g}_n(\mathbf{x})}{\mathrm{\partial }{x}_N}\\ \frac{\mathrm{\partial }{h}_1(\mathbf{x})}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{h}_1(\mathbf{x})}{\mathrm{\partial }{x}_2}& \dots & \frac{\mathrm{\partial }{h}_1(\mathbf{x})}{\mathrm{\partial }{x}_N}\\ ⋮& ⋮& & ⋮\\ \frac{\mathrm{\partial }{h}_m(\mathbf{x})}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{h}_m(\mathbf{x})}{\mathrm{\partial }{x}_2}& \dots & \frac{\mathrm{\partial }{h}_m(\mathbf{x})}{\mathrm{\partial }{x}_N}\end{array}\right].$$

The forward approximation of finite differences provides an estimate of these derivatives, as defined in Eq. (9), where $\mathrm{\Delta }x$ represents the step size and ${\mathbf{e}}_j$ is a unitary vector of the same dimension as x, with a value of 1 for the $j\mathrm{t}\mathrm{h}$ component and 0 for the rest.

(9) $$\mathrm{\nabla }\mathbf{V}(\mathbf{x}{)}_{i,j}\approx \frac{f_i(\mathbf{x}+\mathrm{\Delta }x\cdot {\mathbf{e}}_j)-f_i(\mathbf{x})}{\mathrm{\Delta }x}.$$

This repair method aims to transform x into a feasible solution, which involves setting the elements of the vector $\mathbf{V}(\mathbf{x})$ to zero. By employing Newton-Raphson’s method (Burden, Faires & Burden, 2015), a repaired vector ${\mathbf{x}}_{k+1}$ can be iteratively obtained through Eq. (10), which represents a linear approximation of $\mathbf{V}({\mathbf{x}}_k)$ in the direction of the origin. However, it is common that the number of variables differs from the number of constraints. In this case, the $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k)$ matrix is non-invertible and the Moore-Penrose pseudoinverse (Campbell & Meyer, 2009) must be used, as shown in Eq. (11):

(10) $${\mathbf{x}}_{k+1}={\mathbf{x}}_k-\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k{)}^{-1}\mathbf{V}({\mathbf{x}}_k)$$

(11) $${\mathbf{x}}_{k+1}={\mathbf{x}}_k-\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k{)}^{+}\mathbf{V}({\mathbf{x}}_k)$$where x ${ }_{k+1}$ and x ${ }_k$ represent the updated and current values of the vector x, respectively, $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k{)}^{-1}$ and $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k{)}^{+}$ represents the inverse and the Moore-Penrose pseudoinverse matrix of the gradient matrix $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k)$, respectively. A computationally efficient way of finding $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k{)}^{+}$ is by employing singular value decomposition (Barata & Hussein, 2012). Before proceeding with the pseudoinverse, removing all zero elements of $\mathbf{V}({\mathbf{x}}_k)$ and their corresponding values in $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k)$ is important.

Gradient-based repair method for MINLP problems

As mentioned above, MINLP problems can be analyzed as multiple subproblems, some of which may be highly constrained. This impacts the performance of EAs, whose operators are typically insensitive to constraints, leading to the generation of many infeasible solutions. We propose a basic DE implementation with a gradient-based repair method for MINLP problems, G-DEmi. The repair strategy explores the subproblems independently to improve the exploration inside them. Specifically, for repairing a vector with mixed variables [x, y], only the real variables x are modified while the integer variables y remain fixed.

In contrast to the variants previously described (Chootinan & Chen, 2006; Takahama & Sakai, 2009) with a probability-based approach for repairing solutions, our method repairs only trial vectors u that satisfy two conditions: (i) they lost the tournament against their target vectors, but have a better objective function value, and (ii) they belong to a subproblem (y) where no solution has been repaired in the current generation. These conditions aim to promote the repair of trial vectors with a higher potential to improve OF, exploring each subproblem independently without repairing similar solutions multiple times.

The repair method is explained in Algorithm 1. A mixed trial vector $\mathbf{u}=[{\mathbf{x}}_k^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}}{]}^{\mathrm{T}}$ is defined, where only the ${\mathbf{x}}_k^{\mathrm{u}}$ component is updated during the iterative repair process. As a result, the constraint violation degree vector and the gradient matrix can be found by Eqs. (12) and (13), respectively. Both elements are updated in lines 3 and 4 of Algorithm 1. After removing all zero elements of $\mathbf{V}({\mathbf{x}}_k^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ and their corresponding values in $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$, the pseudoinverse $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_k^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}}{)}^{+}$ is determined in line 6. Finally, the trial vector is updated between lines 7 and 9. The algorithm stops when all elements of vector $\mathbf{V}({\mathbf{x}}_k^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ become zero, indicating a feasible trial vector. Additionally, two stopping criteria are employed to prevent long repair cycles that could cause arbitrary increases in execution time. The first is the maximum number of iterations, which stops the process after $k_{max}$ iterations. The second is the minimum tolerance, which stops the process if the maximum variation in the trial vector during repair is lower than $T_{min}$.

(12) $$\mathbf{V}(\mathbf{u})=\mathbf{V}({\mathbf{x}}_k^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$$

(13) $$\mathrm{\nabla }\mathbf{V}(\mathbf{u})=\frac{({\mathbf{x}}_k^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})}{\mathrm{\partial }{\mathbf{x}}_k^{\mathrm{u}}}.$$

The repairing procedure is illustrated with the following example: consider a scenario with an inequality and an equality constraint, given by Eq. (14).

(14) $$\begin{array}{r}{g}_1(\mathbf{x},y)={x_1}^2+{x_2}^2+y^2-12\le 0\\ h_1(\mathbf{x},y)=x_1+x_2+y-5.5=0.\end{array}$$

Assume $\mathbf{u}=[211{]}^{\mathrm{T}}$ and an equality tolerance of $ϵ=1\times {10}^{-04}$. In the first iteration (where $k=1$), ${\mathbf{x}}_1^{\mathrm{u}}=[21{]}^{\mathrm{T}}$ and ${\mathbf{y}}^{\mathrm{u}}=1$. Therefore, $\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ and $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ are computed by Eqs. (15) and (16), respectively.

(15) $$\mathbf{V}({\mathbf{x}}_1^{\text{u}},{\mathbf{y}}^{\mathbf{\text{u}}})=\left[\begin{array}{c}max(-6,0)\\ -max(1.5-ϵ,0)\end{array}\right]=\left[\begin{array}{c}0\\ -1.4999\end{array}\right]$$

(16) $$\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})=\left[\begin{array}{cc}\frac{\mathrm{\partial }{g}_1({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{g}_1({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})}{\mathrm{\partial }{x}_2}\\ \frac{\mathrm{\partial }{h}_1({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})}{\mathrm{\partial }{x}_1}& \frac{\mathrm{\partial }{h}_1({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})}{\mathrm{\partial }{x}_2}\end{array}\right]=\left[\begin{array}{cc}4& 2\\ 1& 1\end{array}\right]$$

As can be seen, only $h_1({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ was violated. Therefore, the element of $g_1({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ in $\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ and their corresponding values in $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ need to be removed. This leads to $\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})=-1.4999$. Then, $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ and its pseudoinverse $\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}}{)}^{+}$ are computed as shown in Eq. (17).

(17) $$\mathrm{\nabla }\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})=\left[\begin{array}{cc}1& 1\end{array}\right]\Rightarrow \mathrm{\nabla }\mathbf{V}({\mathbf{x}}_1^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}}{)}^{+}=\left[\begin{array}{c}0.5\\ 0.5\end{array}\right]$$

Using Eq. (11), the vector ${\mathbf{x}}_2^{\mathrm{u}}$ is generated as follows:

$${\mathbf{x}}_2^{\mathrm{u}}=\left[\begin{array}{c}2\\ 1\end{array}\right]-\left[\begin{array}{c}0.5\\ 0.5\end{array}\right]\left[\begin{array}{c}-1.4999\end{array}\right]=\left[\begin{array}{c}2.75\\ 1.75\end{array}\right]$$

The updated vector $\mathbf{V}({\mathbf{x}}_2^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ is obtained by Eq. (18). As can be seen, the values of $({\mathbf{x}}_2^{\mathrm{u}},{\mathbf{y}}^{\mathrm{u}})$ satisfy all constraints. Therefore, the trial vector has been successfully repaired, and its new values are $\mathbf{u}=[2.751.751{]}^{\mathrm{T}}$.

(18) $$\mathbf{V}({\mathbf{x}}_2^{\text{u}},{\mathbf{y}}^{\mathbf{u}})=\left[\begin{array}{c}max(-0.3750,0)\\ max(0-ϵ,0)\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

Figure 2 provides an explanation of the gradient-based repair method within the context of DE, for a hypothetical MINLP problem. The feasible region is represented by a gray area, and includes discontinuous feasible parts generated by discrete variables shown as red lines. The objective space is defined by three contour curves: the worst level, a better level, and the best level. In Fig. 2A, the standard DE variant illustrates two instances of target-trial selection. Notably, the target solutions generated infeasible trials, making the target vectors the survivors in both cases, as indicated by the arrow directions. On the other hand, the gradient-based repair DE variant takes into account that both trial solutions have a better value in the objective space than their respective target vectors. This implies that both trial vectors must be repaired and compete again with their targets. Figure 2B illustrates how the repaired vectors enter the feasible regions and outperform their targets. This analysis highlights the advantages of the gradient-based repair method over the traditional selection method for finding feasible solutions within potentially promising subproblems, thereby promoting exploration within better feasible regions. Given the insensitivity of DE to constraints, repair might be essential to achieve feasibility in highly constrained subproblems.

Overall implementation

Algorithm 2 describes G-DEmi. To generate the initial population ${\mathbf{P}}_g$ in step 1, random real and integer values for the solution vector $[\mathbf{x},\mathbf{y}]$ are taken. The objective function $f({\mathbf{x}}_g)$ and the constraint violation degree $G({\mathbf{x}}_g)$ are then evaluated. In step 8, a trial vector ${\mathbf{u}}_{g,i}$ is generated for each target vector ${\mathbf{x}}_{g,i}$, using mutation and binomial crossover (rand/1/bin). The integer variables in ${\mathbf{u}}_{g,i}$ are rounded in step 9 before evaluating the vector in step 10. In the selection operation (steps 11 to 20), the trial vector is compared to its corresponding target vector. In step 11, the best solution is selected according to the feasibility rules described in “Base Methods”. If the trial vector fails to improve its target but still has a lower objective function value $(f({\mathbf{u}}_{g,i})<f({\mathbf{x}}_{g,i}))$, and no other vector of the same subproblem has been repaired in the current generation $({\mathbf{y}}_{g,i}^{\mathrm{u}}\notin \mathbf{Y})$, this trial vector is repaired in step 14. In step 16, following the feasibility rules, the better solution between the repaired vector and its target vector passed to the population of the next generation ${\mathbf{P}}_{g+1}$. Through these steps, G-DEmi generates a new population in each generation.

Settings

Twenty-eight benchmark MINLP minimization problems (F1-F28) were employed to thoroughly evaluate the performance of G-DEmi. They were proposed in Liu et al. (2021a), Datta & Figueira (2013), Liu et al. (2023) and MINLPLib (2023). Some problems have several discontinuous feasible parts with different sizes, and the best solution is located in the smallest part (F1-F4). Other problems have equality constraints (F5, F7, F13, and F14), large dimensionality (F8-F12 and F15, F16, and F18), or binary variables (F17, F22, and F23). Finally, F19 and F26-F28 are highly constrained problems. The details of these test functions can be found in Supplemental File S1, and a summary of their main features is presented in Table 1, where $n$ is the number of variables of the problem; $n_1$ and $n_2$ are the numbers of continuous and discrete variables, respectively; IC and EC are the numbers of inequality and equality constraints, respectively; and AC is the number of active constraints for the best-known solution so far.

A performance evaluation was conducted to compare G-DEmi with six state-of-the-art algorithms. Four different versions of DE were considered, including MI-EDDE (Mohamed et al., 2019), BOToP (Liu et al., 2021b), FROFI (Wang et al., 2015), and SADE (Qin & Suganthan, 2005; Huang, Qin & Suganthan, 2006), as well as two EAs based on different paradigms: PSOmv (Wang, Zhang & Zhou, 2021) based on particle swarm optimization, and EDAmv (Molina-Pérez et al., 2022) based on estimation of distribution. The versions of FROFI and SADE were combined with the cutting and repulsion strategy (CaR) proposed in Liu et al. (2021a), resulting in FROFI-CaR and SADE-CaR, respectively. Furthermore, we conducted experiments to evaluate the effectiveness of the gradient-based repair method in other state-of-the-art DE variants. Two recent approaches were considered for this purpose: chaotic local search-based differential evolution (CJADE) (Gao et al., 2019) and DE with CaR plus surviving solutions (DE-CaR+S) (Molina-Pérez et al., 2024). CJADE has demonstrated high efficacy in optimization problems with continuous variables, while DE-CaR+S has reported remarkable results in MINLP problems. We performed the following comparative analysis: CJADE-CaR (CJADE with CaR strategy) vs. G-CJADEmi (CJADE with gradient-based repair) and DE-CaR+S vs. G-DEmi+S (DE with gradient-based repair plus surviving solutions).

Due to the stochastic nature of EAs, more than simply assessing the quality of outcomes is required to evaluate their performance. Hence, in this experiment, descriptive and inferential statistics were employed. To achieve this, each problem was solved with a maximum of 200,000 OF evaluations and 30 independent runs. The equality-constraint tolerance was set at $1\times {10}^{-04}$. A run was considered successful if $|f({\mathbf{x}}_{best})-f({\mathbf{x}}^{\ast })|\le 1\times {10}^{-04}$, where ${\mathbf{x}}^{\ast }$ denotes the best-known solution and ${\mathbf{x}}_{best}$ represents the best solution generated by the algorithm. All algorithms and their instances were executed using the same computing resources: an Intel Core i7 -4790 CPU @3.6 GHz and 32 GB RAM with Windows 10 Enterprise N and MATLAB 2023a (The MathWorks, Natick, MA, USA).

The iRace parameter tuning tool was used to select parameter sets for all algorithms (López-Ibáñez et al., 2016). The tuning process used a representative set of test problems consisting of F1, F3, F5, F7, F9, F10, F12, F14, F16, F18, F20, F22, F24, F26, and F28. For each algorithm, 2,000 experiments were conducted in iRace. The parameters used in each algorithm are presented in Table 2, where population size is denoted by NP, crossover rate by CR, scaling factor by F, maximum number of repair iterations by $k_{\mathrm{m}\mathrm{a}\mathrm{x}}$, minimum repair tolerance by $T_{\mathrm{m}\mathrm{i}\mathrm{n}}$, proportion of population partitions by $P_a$, and initial tolerance for equality constraints by $e_0$. Parameters $t_c$ and $c_p$ are the controllers of the e-constrained method. For algorithms with a learning process, L denotes the learning iteration count, $\alpha$ is the learning rate, $r$ represents the chaotic search radius, $c$ is the rate of parameter adaptation, $p$ is the greediness of the mutation strategy, $ϵ$ is the link parameter, $r_M$ is the mutation rate, $\beta$ represents the mutation factor, $w$ denotes the number of histogram bins, $ϵ$ represents the end bins parameter, $A_c$ is the acceleration coefficient, and $\rho$ represents the shrinking parameter. The failure threshold is denoted by T. In replacement strategies, NR indicates the number of surviving replacements, and $\delta$ is the improvement tolerance. Additionally, $l_1$ and $l_2$ represent the count of weighted difference vectors.

As described above, G-DEmi rounds the values of integer variables after the rand/1/bin operation. The value of F suggested by iRace (higher than the 0.5 recommended by Storn & Price, 1997) can be effective in reducing the repeated integer values after rounding and encouraging the exploration of different discontinuous feasible parts. The value of CR indicates that crossover is also significant in this scheme to prevent accelerated convergence. The population of 50 individuals coincides with several successful DE experiments for MINLP problems (Mohamed et al., 2019; Liu et al., 2021a).

Benchmark results and analysis

The experimental results were analyzed using descriptive statistics. The selected indicators are widely recognized in the field for assessing algorithms, providing a comprehensive view of their performance (Liu et al., 2021b, 2021a). The feasible rate (FR) represents the percentage of runs in which the algorithm finds at least one feasible solution, whereas the successful rate (SR) represents the percentage of runs in which the algorithm finds the optimal solution within the given tolerance. Ave and Std denote the average and the standard deviation of the best OF values of each problem from its 30 independent runs. If an algorithm cannot achieve 100 $\mathrm{\%}$ FR, then the Ave and Std values are denoted as “NA”. Mean FR and mean SR denote the mean feasible rate and the mean successful rate achieved by the algorithm over the twenty-eight test problems. T represents the average execution time for each problem in seconds.

Inferential statistics were employed to evaluate the statistical significance of the differences in results. A comparison between G-DEmi and the other algorithms was carried out using the Wilcoxon Rank-Sum Test (WRST) with a 0.05 significance level (Derrac et al., 2011). Additionally, we applied the Bonferroni correction to control the family-wise error rate that arises from multiple pairwise comparisons (Dinno, 2015). The symbols indicating the superiority, inferiority, or similarity of our proposal compared to each competitor are denoted by “+,” “−,” and “ $\approx$”, respectively.

Tables 3 and 4 present the statistical values of each problem, in separate rows. For example, in problem F1, all algorithms successfully found feasible solutions across all executions, yielding 100% FR. The success rates indicate that MI-EDDE failed every run (SR = 0), whereas the other competing algorithms consistently achieved the optimal solution (SR = 100). Based on WRST, G-DEmi was comparable to the other competitors, except for MI-EDDE, which was significantly inferior. Additionally, MI-EDDE had the fastest average execution time of 5.91 s. The remaining problems were analyzed similarly.

Mean FR, mean SR, and the final count of WRST are provided in the last section of Table 4. It can be observed that G-DEmi reached feasible solutions in all executions for all problems, 100% of Mean FR, surpassing the 84.52% for MI-EDDE, 82.86% for BOToP, 85.83% for SADE-CaR, and 90.36% for FROFI-CaR. Note that the highly constrained condition of problems F19 and F26-F28 had no influence on the FR of G-DEmi. However, this condition significantly affected the performance of the competing algorithms. In terms of mean SR, G-DEmi obtained the highest value at 98.10%, followed by FROFI-CaR at 73.45%, SADE-CaR at 59.52%, BOToP at 46.79%, and MI-EDDE at 38.10%. The final count of WRST shows that G-DEmi significantly outperformed MI-EDDE in 20 problems, whereas both algorithms had similar results in eight problems, with MI-EDDE never surpassing the performance of G-DEmi. BOToP is significantly outperformed by G-DEmi in 17 problems, equaled in 11 problems, and never achieving better results than G-DEmi. Similarly, FROFI-CaR was surpassed by G-DEmi in 10 problems, whereas its results were similar for 17 problems, and it only significantly exceeded G-DEmi in one problem (F11).

In this experiment, there were some outlier values that should be noted. For problems F1-F3, the solutions of MI-EDDE tended to converge towards local optima deliberately located in the most accessible discontinuous feasible parts. In contrast, problems F7 and F13, with equality constraints, affected the performance of other algorithms but not the G-DEmi proposal. In problems F8-F10, certain points were attractive to the algorithm population, with MI-EDDE and POToP showing the worst performances. Problem F11 has more than 10 real variables and a high sensitivity of the objective function to small variations in the decision variables. As a result, G-DEmi did not reach the solution before finishing the execution in 40% of the cases. Remarkable advantages of G-DEmi over PSOmv and EDAmv are shown in Tables S5 and S-6 in Supplementary File S2. PSOmv and EDAmv obtained Mean FR values of 85.83% and 88.45%, and Mean SR values of 14.17% and 37.74%, respectively. G-DEmi outperformed PSOmv in 25 problems, with only three instances of similar performance. On the other hand, EDAmv was outperformed in 21 problems and only achieved similar performances in seven problems. G-DEmi did not exhibit inferior performance with respect to PSOmv and EDAmv in any case.

The experiments involving the gradient-based repair method with other DE variants are reported in Tables 5 and 6. These tables present the results of CJADE-CaR vs. G-CJADEmi and DE-CaR+S vs. G-DEmi+S. The comparison highlights that G-CJADEmi outperformed CJADE-CaR, achieving 100% of Mean FR and 90.48% of Mean SR compared to 98.81% and 69.05%, respectively. G-CJADEmi exhibited better results in 12 problems and similar results in 14 problems, but CJADE-CaR performed better in problems F15 and F23. In the second competition, G-DEmi+S surpassed DE-CaR+S with a Mean FR of 99.52% vs. 90.24% and a Mean SR of 91.55% vs. 83.93%. G-DEmi+S achieved better results in seven problems, with no significant difference in 19 problems, and was outperformed in F12 and F23 by DE-CaR+S. In general, the variants with gradient-based repair exhibited notable superiority, with G-DEmi as the most promising in terms of feasibility and success.

Additionally, an analysis was conducted to examine the convergence of G-DEmi during the evolutionary process. Three problems were selected for this: F3 with two variables, F12 with 15 variables and high sensitivity of OF to changes in decision variables, and F27 highly constrained. Figure 3 shows the convergence curves of the median solutions in terms of solution quality for problem F12. Each curve starts from the generation of at least one feasible solution. A star marker on each curve shows when the algorithm reached the global solution. It is observed that G-DEmi converged rapidly and found the global optimum before 100,000 evaluations (blue point). In contrast, the other algorithms consumed 200,000 evaluations without achieving this result. In Supplemental File S2, Figs. S1 and S2 show the curves for F3 and F27, indicating an even more favorable outcome for G-DEmi.

Likewise, the impact of the parameters $k_{max}$ and $T_{min}$ on the performance of G-DEmi was analyzed through the described problems. For each problem, 30 executions of 200,000 evaluations were carried out using different combinations of them, while the other parameters remained fixed at previously reported values. In Supplemental File S2, Tables S2, S3, and S4 show the average execution times [s] and the average OF values for each instance. The best values are highlighted in bold. In general, the performance was adversely affected by the option of no repair ( $k_{max}=0$), and better performances were seen as $k_{max}$ increased and $T_{min}$ decreased. However, it is important to note that a longer repair process resulted in higher time consumption. Therefore, the execution times could be a determining factor for reasonable values such as those used. Values of $k_{max}$ equal to 50 and 100, and $T_{min}$ of $1\times {10}^{-60}$ and $1\times {10}^{-80}$ proved effective in this study. The consistency of the performance suggests that both parameters have a robust behavior concerning the algorithm’s performance.