R/constraint_vbr.R

In MOEADr: Component-Wise MOEA/D Implementation

#' "Violation-based Ranking" constraint handling method for MOEA/D
#'
#' Uses the Violation-based Ranking handling method to generate a
#' preference index for the MOEADr framework.
#'
#' This function calculates the preference index of a set of neighborhoods
#' based on the "violation-based ranking" (VBR) constraint handling method. Please
#' see [order_neighborhood()] for more information on the preference index
#' matrix.
#'
#' The VBR strategy generalizes some well-known methods for handling constraints
#' in population-based metaheuristics (see Section `c(x) Criteria`).
#' This strategy essentially ranks points within for a given subproblem based on
#' their aggregated function value (`f^{agg}(x|w_i)`) or their total constraint
#' violation (`v(x)`). Specific variations of this strategy differ on the
#' criteria for using one or the other.
#'
#' The value used for ranking a given point `x` can be summarized as:
#'
#'\tabular{lll}{
#'   Violation  \tab | c(x) criterion \tab | Rank using: \cr
#'   `v(x) = 0` \tab | `c(x) = *`     \tab | `f^{agg}(x|w_i)`\cr
#'   `v(x) > 0` \tab | `c(x) == TRUE` \tab | `f^{agg}(x|w_i)`\cr
#'   `v(x) > 0` \tab | `c(x) == FALSE`\tab | `v(x)`      \cr
#'}
#'
#' Points compared according to their `f^{agg}(x|w_i)` values (i.e., feasible
#' points and those for which `c(x) = TRUE`) are ranked first (i.e., receive
#' ranks between `1` and `n_{feas}`, where `n_{feas}` is the
#' number of feasible points in the i-th neighborhood), with points that are
#' compared according to their `v(x)` values receiving ranks between
#' `(n_{feas} + 1)` and `T + 1` (`T` being the size of the neighborhood. The `+1`
#' comes from including the incumbent solution in the comparison).
#'
#' @section c(x) Criteria:
#' Specific variations of the VBR differ on how the criterion c(x) is
#' implemented. Three variants are currently implemented in the MOEADr package:
#'
#'\tabular{lll}{
#'   Method                                 \tab | ID            \tab | `c(x)` \cr
#'   Tournament Selection `[Deb2000]`       \tab | `$type = "ts"`\tab | `FALSE`\cr
#'   Stochastic Ranking `[Runarsson2000]`   \tab | `$type = "sr"`\tab | `runif() < pf`\cr
#'   Violation Threshold `[Asafuddoula2014]`\tab | `$type = "vt"`\tab | `v(x) < eps_v^i`\cr
#'}
#'
#' where \eqn{pf \in [0,1]} is a user-defined parameter for the "sr" method, and
#' `eps_v^i` is subproblem-dependent, adaptive quantity calculated internally
#' in the routine (see `[Asafuddoula2014]` and `[Campelo2017]` for details).
#'
#' @section Using an External Archive:
#' For types "sr" and "vt", it is possible for the algorithm to lose feasible
#' solutions during its update step, since there is a non-zero probability of
#' unfeasible solutions replacing feasible ones. In these cases, it is
#' recommended to set the [moead()] parameter `update$UseArchive = TRUE`, so
#' that an external archive is built with the best feasible solutions found for
#' each subproblem.
#'
#' @param bigZ Matrix of scalarized objective values for each neighborhood and
#' the incumbent solution (generated by [scalarize_values()])
#' @param bigV Matrix of violation values for each neighborhood and the
#' incumbent solution (generated in [order_neighborhood()])
#' @param type type of `c(x)` function to use (see `c(x) Criteria` for details).
#' @param pf probability parameter for type = "sr" (ignored in other modes).
#' @param ... other parameters (unused, included for compatibility with
#' generic call)
#'
#' @return `[ N x (T+1) ]` matrix of preference indices. Each row `i` contains
#' a permutation of `{1, 2, ..., (T+1)}`, where `1,...,T` correspond
#' to the solutions contained in the neighborhood of the i-th subproblem,
#' `B[i, ]`, and `T+1` corresponds to the incumbent solution for that
#' subproblem. The order of the permutation is defined by the specific strategy
#' defined by the input variable `type`).
#'
#' @export
#'
#' @section References:
#' `[Deb2000]` K. Deb,
#' "An efficient constraint handling method for genetic algorithm",
#' Computer Methods in Applied Mechanics and Engineering 186(2–4):311–338, 2000.
#'
#' `[Runarsson2000]` T. Runarsson, X. Yao,
#' "Stochastic ranking for constrained evolutionary optimization",
#' IEEE Transactions on Evolutionary Computation4(3):284–294, 2000.\cr
#'
#' `[Asafuddoula2014]` M. Asafuddoula, T. Ray, R. Sarker, K. Alam,
#' "An adaptive constraint handling approach embedded MOEA/D,”
#' 2012 IEEE Congress on Evolutionary Computation (CEC).\cr
#'
#' `[Campelo2017]`  F. Campelo, L.S. Batista, C. Aranha (2020): The {MOEADr}
#' Package: A Component-Based Framework for Multiobjective Evolutionary
#' Algorithms Based on Decomposition. Journal of Statistical Software
#' \doi{10.18637/jss.v092.i06}\cr
#'
#'

constraint_vbr <- function(bigZ, bigV, type = c("ts", "sr", "vt"), pf = NULL, ...)
{
  # ========== Error catching and default value definitions
  valid_types <- c("ts", "sr", "vt")
  assertthat::assert_that(
    identical(dim(bigZ), dim(bigV)),
    type %in% valid_types
    )

  if (type == "sr"){
    assertthat::assert_that(is.numeric(pf),
                            is_within(pf, 0, 1, strict = FALSE))
  }
  # ==========

  # Feasible points in each neighborhood
  feasible <- (bigV == 0)
  cx <- switch(type,
        # TS: c(x) = FALSE \forall x
        ts = (bigV > Inf),
        # SR: c(x) = (runif() <= pf)
        sr = (randM(bigV) <= pf),
        # VT: c(x) = (v(x) <= eps_v)
        vt = (bigV <= (colSums(feasible) / (nrow(bigV) ^ 2)) * colSums(bigV)))

  # Points to be compared using f^{agg}
  useF <- feasible | cx

  # Create the matrix of performance for feasible indexes,
  # and of violation for infeasible indexes.
  bigF        <- bigZ
  bigI        <- bigV
  bigF[!useF] <- NA
  bigI[useF]  <- NA

  # Sort the feasible and infeasible matrixes, putting
  # all NAs in the back or front, respectively
  indxF <- t(apply(bigF,
                   MARGIN  = 2,
                   FUN     = order,
                   na.last = TRUE))

  indxI <- t(apply(bigI,
                    MARGIN  = 2,
                    FUN     = order,
                    na.last = FALSE))

  # Merge feasible and unfeasible matrixes
  indx.joint <- t(sapply(1:ncol(bigZ),
                         function(i) {
                           (indxI[i, ] * !is.na(bigI[indxI[i, ], i])) +
                             (indxF[i,]   * !is.na(bigF[indxF[i, ], i])) }
                         ))

  return(indx.joint)
}