Routine for real variable SBX crossover taken from the author's implementation (Kalyanmoy Deb)

sbx.c

/* Routine for real variable SBX crossover taken from the author's implementation (Kalyanmoy Deb).
 *
 * Paper describing the algorithm: 
 * Title: An Efficient Constraint Handling Method for Genetic Algorithms
 * Author: Kalyanmoy Deb
 * More info: Appendix A. Page 30.
 * URL: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.33.7291&rep=rep1&type=pdf
 *
 * Commentaries by Tiago Peres França.
 *
 * OBS1: The paper describes the SBX algorithm for boundless variables and variables in a well determined interval. 
 * The following code implements the version with limited bounds.
 *
 * OBS2: Note that the original code wasn't changed and because of the font used here, y1 and yl is a bad choice of
 * variable names. Do not mistake one for the other.
 *
 * OBS3: The implementation here is slightly different from what is described in the paper. In the paper, for 
 * calculating the children, it is used only one value of beta: beta = 1 + (2 / (y2 - y1)) * min(x1 - xl, xu - x2).
 * Here, it is used two values for beta: one for calculating the 1st child and another for calculating the 2nd. The 1st
 * value of beta is calculated as: beta = 1 + (2 / (y2 - y1)) * (min(x1, x2) - xl). The second value of beta is 
 * calculated as: beta = 1 + (2 / (y2 - y1)) * (xu - max(x1, x2)). When y2 - y1 = 0, the division can not be made and
 * the children get the same value as the parents.
 *
 * Constants, functions and global variables definitions (ordered as they appear in the code):
 * - randomperc(): fetches a random number between 0.0 and 1.0;
 * - pcross_real: crossover rate (global);
 * - nrealcross: counts the number of crossovers (global);
 * - nreal: total number of variables in each individual/solution (global);
 * - EPS: precision error tolerance, its value is 1.0e-14 (global constant);
 * - min_realvar: vector of size nreal containing the minimum value for each of the variables in a solution;
 * - max_realvar: vector of size nreal containing the maximum value for each of the variables in a solution;
 * - eta_c: represents "nc" in the paper, ie. the distribution index for the crossover. In the original NSGA-II implementation, eta_c is given by the user.
 */

void realcross (individual *parent1, individual *parent2, individual *child1, individual *child2)
{
    int i;
    double rand;
    double y1, y2, yl, yu; // y1 stores the value for the 1st child; y2 the value for the 2nd child; yl (notice it's an L not a 1) holds the lower limit for the variable; yu holds the upper limit.
    double alpha, beta, betaq; // betaq, in the paper, is the symbol beta with a line above
    if (randomperc() <= pcross_real) // roll the dices. Should we apply a crossover?
    {
        nrealcross++; // increments the number of crossovers
        for (i=0; i<nreal; i++) // for each variable of a solution (individual)
        {
            if (randomperc() <= 0.5 ) // according to the paper, each variable in a solution has a 50% chance of changing its value. This should be removed when dealing with one-dimensional solutions.
            {
                // the following if/else block puts the lowest value between parent1 and parent2 in y1 and the other in y2.
                if (parent1->xreal[i] < parent2->xreal[i])
                {
                    y1 = parent1->xreal[i];
                    y2 = parent2->xreal[i];
                }
                else
                {
                    y1 = parent2->xreal[i];
                    y2 = parent1->xreal[i];
                }
                if (fabs(parent1->xreal[i]-parent2->xreal[i]) > EPS) // if the value in parent1 is not the same of parent2
                {
                    yl = min_realvar[i]; // lower limit of the i-th variable of an individual. Note that yl != y1.
                    yu = max_realvar[i]; // upper limit of the i-th variable of an individual
                    // CALCULATION OF THE FIRST CHILD
                    rand = randomperc();
                    beta = 1.0 + (2.0*(y1-yl)/(y2-y1)); // it differs from the paper here. The paper uses one value of beta for calculating both children. Here, we use one beta for each child.
                    alpha = 2.0 - pow(beta,-(eta_c+1.0)); // calculation of alpha as described in the paper
                    // calculation of betaq as described in the paper
                    if (rand <= (1.0/alpha))
                    {
                        betaq = pow ((rand*alpha),(1.0/(eta_c+1.0)));
                    }
                    else
                    {
                        betaq = pow ((1.0/(2.0 - rand*alpha)),(1.0/(eta_c+1.0)));
                    }
                    child1->xreal[i] = 0.5*((y1+y2)-betaq*(y2-y1)); // calculation of the first child as described in the paper
                    // CALCULATION OF THE SECOND CHILD
                    beta = 1.0 + (2.0*(yu-y2)/(y2-y1)); // differs from the paper. The second value of beta uses the upper limit (yu) and the maximum between parent1 and parent2 (y2)
                    alpha = 2.0 - pow(beta,-(eta_c+1.0)); // calculation of alpha as described in the paper
                    // calculation of betaq as described in the paper
                    if (rand <= (1.0/alpha))
                    {
                        betaq = pow ((rand*alpha),(1.0/(eta_c+1.0)));
                    }
                    else
                    {
                        betaq = pow ((1.0/(2.0 - rand*alpha)),(1.0/(eta_c+1.0)));
                    }
                    child2->xreal[i] = 0.5*((y1+y2)+betaq*(y2-y1)); // calculation of the second child as described in the paper
                    // ensures the value of both children are in the correct bounds [yl, yu]. According to the paper, this should not be needed.
                    if (child1->xreal[i]<yl)
                    {
                        child1->xreal[i]=yl;
                    }
                    if (child1->xreal[i]>yu)
                    {
                        child1->xreal[i]=yu;
                    }
                    if (child2->xreal[i]<yl)
                    {
                        child2->xreal[i]=yl;
                    }
                    if (child2->xreal[i]>yu)
                    {
                        child2->xreal[i]=yu;
                    }
                }
                // the paper is not very clear about this, but I assume, in the equation of beta (not betaq), y2 and y1, since they could not have been calculated yet, refer to the parents. So, if both parents are equal at the specified variable, the divisor would be zero. In this case, the children should have the same value as the parents. 
                else // if the i-th variable has the same value in both parents
                {
                    child1->xreal[i] = parent1->xreal[i];
                    child2->xreal[i] = parent2->xreal[i];
                }
            }
            else // 50% chance of changing values. In the case random > 0.5, the children should have the same value as the parents
            {
                child1->xreal[i] = parent1->xreal[i];
                child2->xreal[i] = parent2->xreal[i];
            }
        }
    }
    // if the random number generated is greater than the crossover rate, return the children as exact clones of the parents
    else
    {
        for (i=0; i<nreal; i++)
        {
            child1->xreal[i] = parent1->xreal[i];
            child2->xreal[i] = parent2->xreal[i];
        }
    }
    return;
}

Based on the following paper: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.33.7291&rep=rep1&type=pdf
Functional implementation in Javascript: https://gist.github.com/Tiagoperes/6849da791bf189bb2746e010b87b2c2f

Fala Thiago, meu projeto de fim de curso envolve o uso de algoritmo genético e o artigo no qual me baseio utiliza desse método de cruzamento e não existe tanto material online disponível sobre, e o seu código está me ajudando a elucidar dúvidas que eu vinha tendo sobre este método. E UFU ao contrário é UFU!!

Haha. Vlw mtfrsantos!

Thank you! This was super helpful!

@Tiagoperes Sua implementação esta melhor que a do artigo.
Eu estava usando essa implementação aqui https://rdrr.io/cran/nsga2R/man/boundedSBXover.html, que como vc disse, utiliza um beta só com o min((y1 - yL), (yu - y2)). O problema desta implementação é que os filhos vão ter a mesma "folga" mínima entre o valor dos pais e os extremos.
Por exemplo: Se os pais tem valores 5 e 15 e os limites são 0 e 100, a folga para o mínimo é 5 e a folga pro máximo é 85, utilizando a versão com apenas um beta, a folga tanto para o máximo quanto para o mínimo vai ser a mínima entre as duas folgas, ou seja, 5.
Estas figuras ilustram o problema.
Versão com um beta utilizando min((y1 - yL), (yu - y2)):

Versão com dois betas, beta1 = 1.0 + (2.0*(yu-y2)/(y2-y1)); e beta2 = 1.0 + (2.0*(y1-yl)/(y2-y1));

Nestes exemplos foi removida a probabilidade de 50% de não alterar a variável.

Obrigado @Tiagoperes.