2.1 Adaptive-width Histogram Model

The AWH model promotes promising regions by assigning them high probabilities, while in the other regions very low probabilities are assigned. One AWH is developed for each decision variable independently.

The search space [ai,bi] of the ith variable xi is divided into (W+2) bins (regions), to define the probabilities Pric for the AWH model.

Points [pi,0,pi,1,…,pi,w+1,pi,w+2] define the width of the bins shown in Fig. 2, where pi,0=ai and pi,w+2=bi (ai and bi are the lower and upper bounds of xi, respectively).

Fig. 2 Search progress of the AWH model for W = 3. (a) first generations, (b) later generations 

The total number of bins is (W+2) although the input parameter for EDAmv is W , since the algorithm creates two more bins: one between the lower boundary ai and the point pi,1, and another one between the point pi,w+1 and the upper boundary bi (unpromising regions).

By assuming that xi,min⁡1 and xi,min⁡2 are the smallest and the second smallest existing values of variable xi, respectively, and xi,max⁡1 and xi,max⁡2 are the highest and the second highest existing values of variable xi, respectively, then points pi,1 and pi,w+1 are defined as in (7) and (8):

The W bins of the promising areas are located in the range [pi,1,pi,w+1], and have the same width a, given by (9):

Let Ai,j be the count of individuals for the ith variable located in the jth bin. As can be seen in Fig. 2, the end bins do not contain solutions (unpromising regions), then Ai,1=Ai,W+2=.

However, a small value will be assigned through the parameter eb, to avoid premature convergence. Ai,j is obtained by (10):

The first case in (10) is the count of bins of promising regions [pi,1,pi,w+1]. The second case corresponds to unpromising regions with ebvalue.

The third case assigns zero to the end bins with empty range. The probability of the ith variable in the jth bin is obtained by (11):

2.2 Learning-based Histogram Model Linked with ε-constrained

The LBHε model is used for handling integer variables. It is a link between the LBH model and the ε-constrained method.

The aim is to maintain an equal probability for all available integer values until ε reaches a predefined value εp, as is shown in Fig. 3 (a).

Fig. 3 LBHε model for v = 6. (a) ε>εp equal probability for all available integer values, (b) ε≤εp considering population distribution 

When ε reaches εp, the LBH model begins the learning process, i.e., considering the information of the population distribution to update the probability, as shown in Fig. 3 (b).

If the ε-constrained method has been effective, for values of ε sufficiently small, the solutions must be close to those parts of the feasible region with promising objective function values.

Therefore, if the histogram begins the learning process at that point, it has a better chance of converging to good solutions.

Considering that the variable ym has v available integer values, with v∈{Lm,Lm+1,Lm+2,…,Um}, the probability of the nth available value of v is defined by (12):

where N is the population size, t is the current generation, γ is the population learning rate, and Countv is the number of individuals with the nth available value of v.

Let lmax be the maximum number of generations, and γ a dynamic parameter defined by (13):

Therefore, as the number of generations advances, γ gradually increases as well, which implies an accelerated learning process, i.e., the model uses more information of the current population distribution.

2.3 Sampling

After the histograms have been developed, the offspring is obtained by sampling the models.

In case of a continuous variable xi, a bin j is firstly selected according to a randomly generated probability, then xi is uniformly sampled from the points that limit the bin selected [pi,j−1,pi,j).

For a discrete variable ym, an available value of v∈{Lm,Um} is selected by a randomly generated probability.

2.4 Hybridization with a Mutation Operator

The mutation operation is added to generate the real variables. The vector of real variables x of each offspring is generated by mutation or by sampling taking into account the predefined mutation probability rM, i.e. if this probability is satisfied for a solution vector, its real variables are computed as shown in (14) and (15):

where k,i are the index of the current solution vector and current variable, respectively, g is the current generation, xbest,ig is the ith variable of the best solution vector found so far, randk,i is a random number between 0 and 1, and βmin⁡ and βmax⁡ are the lower and upper bounds of β predefined by the user, with values between 0 and 1.

In the new proposal the values of βmin⁡ and βmax⁡ will always be set to 0 and 1, respectively.

2.5 Constraint Handling

The replacement mechanism to get the next population is carried out through parent-offspring competition using the ε-constrained method.

The ε-constrained method was proposed by Takahama and Sakai [¹⁴] as a constraint-handling technique.

Given two function values f(x1), f(x2), and two constraint violations ϕ(x1), ϕ(x2) for two points x1 and x2, the ε-constrained method uses the ε-level comparisons described in (16):

where ε-level comparisons are defined as an order relation on a pair of objective function and constraint violation values (f(x),ϕ(x)).

This means that the candidates with a violation sum lower than ε are considered as feasible solutions and are ordered according to their fitness values.

In the case of ε=0, ϕ(x) always precedes f(x). Therefore, this method favors the approach to the feasible region by keeping slightly infeasible solutions with promising fitness values.

The ε-level decreases at each iteration G until the predefined iteration number Tc is reached, after that ε=0, as indicated by (18):

where cp is a parameter to control the speed of constraint relaxation, ε(0) is the initial value of ε, and xθ is the top θth in an array sorted by total constraint violation (θ=0.2N).

3.1 LBHε Improvement

As described in (12), γ is the learning rate of the population. A high value of γ increases the role of the population distribution Countv/N in obtaining the Prmd(t), whereas a low value mainly considers the histogram of the previous generation, Prmd(t−1).

However, when certain admissible values begin to prevail statistically over others, the histograms and the populations begin to be similar, so the terms of the equation (12), instead of combining different information, emphasize the same search direction and cause accelerated (and often premature) convergence.

In this work, the following LBHε model is proposed:

where Pemd are equal probabilities for all v values of the mth variable, and are given by (20):

In this model, Pemd contributes to the exploration of the algorithm, while Countv/N contributes to the exploitation on the most populated regions (promising regions).

As can be seen in Fig. 4, now for very low values of γ, the histogram will be flatter (low selection pressure).

Fig. 4 LBHε model, γ=0 random exploration, γ=0.5 middle consideration of population distribution, γ=1 total consideration of population distribution 

As the value of γ increases, the histogram and selection pressure will be more consistent with the population distribution. As in the previous case, γ is a dynamic parameter defined by (13).

3.2 Repulsion

The repulsion strategy proposed in [⁶] consists of two steps: (i) judge whether the population is trapped into a solution, and (ii) apply a repulsion operator to the discontinuous feasible part containing the solution, and restart the population. Eq. (21) is the fail consideration to find a better solution:

where fbest and gbest are the objective function value and the degree of constraint violation of the best solution found so far, respectively, f′best and g′best are the objective function value and the degree of constraint violation of the best solution in the current generation, respectively.

If (21) is satisfied, it means that the algorithm fails to find a better solution, then the counter is incremented (ctr=ctr+1). If (21) is not satisfied in any generation, the counter is reset (ctr=0).

If ctr is greater than a predefined failure threshold T, the population is considered to be trapped in a solution, and the discontinuous feasible part (y) containing that solution has been explored. Then the population is regenerated, and the solution is recorded in the store archive. Any population member that has a vector y contained in store will be penalized with an arbitrarily large degree of constraint violation.

ε-constrained method is also restarted but with a new Tc value with fewer generations, called fast generation control (T′c). At the end of the execution, the recorded solutions should be considered to return the best solution.

4.1 Benchmark Problems

Sixteen MINLP problems (F1-F16) were used to evaluate the performance of EDAIImv. Because of the space limitation, a detailed description of the problems is not included, but it can be found in [⁶].

The maximum number of objective function evaluations was set at 200,000, and 25 independent runs were executed for each problem. The tolerance value for the equality constraints was set at 1.0E-04.

A run was considered as successful if: |f(xbest)−f(x*)| ≤1.0E-4, where x* is the best known solution and xbest is the best solution provided by the algorithm.

4.2 Algorithms and Parameter Settings

PSOmv [¹⁶], EDAmv [⁹], and EDAIImv were the competing algorithms in the experiment. PSOmv also uses the LBH model for handling discrete variables. However, the γ is an adaptive parameter, and the LBH probability is updated using only the best half of the swarm.

To prove the individual contribution of each modification proposed, the instance with only LBHε improvement (EDAmv(I)) was also included. The algorithms were tuned using the iRace parameter tuning tool [⁸]. The parameter values were as follows:

– PSOmv: swarm size N=300, acceleration coefficient c=1.5299, learning rate γ=0.0125.

– EDAmv: N=50, numbers of bins W=4, end bins parameter eb=2.3959, control generation Tc=3,000, control speed parameter cp=8, link parameter εp=0.2399, and mutation parameters: rM=0.6, βmin⁡=0.3, βmax⁡=0.9.

– EDAmv(I): N=50, W=3, eb=2, Tc=2000, cp=7, εp=5, and rM=0.3.

– EDAIImv: N=50, W=3, eb=2, Tc=2000, cp=7, εp=5, rM=0.3, failure threshold T=400, and fast control generation T′c=200.

4.3 Analysis of Results

Table 1 summarizes the results of PSOmv, EDAmv, EDAmv(I), and EDAIImv.

These results are assessed considering the terms Feasible Rate (FR), Successful Rate (SR), Average (Ave), and Standard Deviation (Std Dev), over 25 independent runs. “NA” means that an algorithm cannot achieve 100% FR.

EDAmv(I) beats EDAmv in nine problems (F7:F12, F14:F16) in at least one of the term concerned, proving that LBHε has a positive influence on the algorithm performance.

As mentioned above, the LBH model (used in EDAmv) has two terms that could contain redundant information, producing an accelerated convergence.

However, for problems F1, F3, and F5, where the solutions are in small feasible parts, a slower convergence of LBHε (used in EDAmv(I)) causes that the ε-level reaches zero value when the histogram has not yet converged to the small promising part.

The repulsion strategy is very useful for this situation, since restarts the exploration in the remaining unexplored parts.

As can be seen, the implementation of repulsion strategy in EDAIImv improves the performance for problems F1, F3, and F5 without compromising the rest of the problems.

A Wilcoxon’s rank-sum test at a 0.05 significance level was carried out between EDAIImv and each competitor, in order to evaluate the significant differences in the results.

In Table 1, [+], [≈] and [−] denote that EDAIImv is better than, worse than, and similar to its current competitor, respectively.

As shown in the final part of Table 1, the results of EDAIImv are significantly better than EDAmv in eight problems (F7:F12, F15,F16), similar in another eight problems (F1:F6, F13, F14), and in no case EDAmv surpasses the result of the new proposal.

EDAmv(I) results are outperformed on five problems (F1, F3, F5, F7, F10), matched on eleven problems (F2, F4, F6, F8, F9, F11:F16), and in no case is EDAIImv outperformed by EDAmv(I).

It is clear that EDAIImv has significantly better results than previous variants. Analyzing the results of this sequenced implementation, it can be concluded that each proposed modification contributes to a better performance.

Regarding PSOmv, the new proposal is significantly better in seven test problems (F1, F3, F5, F7, F11, F12, F14) and no difference in four problems (F2, F4, F6, F13), while PSOmv outperformed EDAIImv in five problems (F8, F9, F10, F15, F16).

Although in general EDAIImv has a better performance than PSOmv, the advantage of PSOmv in the last mentioned problems is due to a superior diversity in the exploration.

Therefore, it is recommended in future works to focus on promoting greater diversity in EDAIImv.