R: Why does the optimization solution from mco's nsga2 function not satisfy the constraints?

Question

Below is a reproducible example. First my data:

myvec <- c(-0.531346314931298, 0.298659434001583, -0.691720475359108)
mymat <- new("dgeMatrix", x = c(16.3574344575224, -15.8734447380907, -2.5840252937711, 
                                3.06801501320276, 1.02979983522784, -0.54581011579616, -18.9562379944301, 
                                -13.26737180927, -15.3314450728231, -16.892164730877, -16.232802331708, 
                                -15.9908074719921, 8.05590245092502, -8.05590245092502, -1.51374510440886, 
                                1.51374510440887, 0.12099742985792, -0.120997429857927, -9.56964755533389, 
                                -6.54215734651616, -8.05590245092502, -8.05590245092502, -7.9349050210671, 
                                -8.17689988078295, 3.78595636574668, -4.26994608517835, -1.09039597003911, 
                                0.606406250607442, -0.632174774229305, 0.14818505479763, 1.33542428436906, 
                                -0.303932077216888, 4.33320538312598, -3.30171317597381, 25.7164288780753, 
                                -24.6849366709232, 2.96628530035713, -4.59014893387888, 24.3728836566943, 
                                -25.9967472902161, -1.06980486854892, -0.554058764972834, 25.5514562424382, 
                                24.7395244256774, 28.9907656065072, 21.3002150616084, 24.8998122739105, 
                                25.3911683942051, 2.74913208137828, -1.12526844785653, 13.3846769837898, 
                                -11.760813350268, 1.33987291197775, 0.283990721544011, -0.921712011956532, 
                                -0.109780195195655, 1.3543960057099, -2.38588821286208, 12.0963242028502, 
                                -13.1278164100024, -12.8382805804706, -12.4583324395663, -14.9147553487727, 
                                -10.3818576712642, -13.0374407869434, -12.2591722330935, -1.30406270512594, 
                                1.30406270512594, 1.49880108324396e-15, -1.77635683940025e-15, 
                                6.32415325500923, -6.32415325500923, -1.11022302462516e-15, 2.33146835171283e-15, 
                                -1.30406270512594, 1.30406270512594, 6.32415325500923, -6.32415325500922, 
                                -1.30406270512594, 1.30406270512594, -5.55111512312578e-16, 0, 
                                6.32415325500923, -6.32415325500923, -6.32415325500922, -6.32415325500923, 
                                -7.62821596013516, -5.02009054988329, -6.32415325500923, -6.32415325500923
), Dim = c(30L, 3L), Dimnames = list(NULL, NULL), factors = list())
myvec2 <- new("dgeMatrix", x = c(2.43138842392909e-14, 1.77635683940025e-14, 
                                 -0.0990891178849278, 0.0990891178849509, 1.4210854715202e-14, 
                                 -1.99840144432528e-15, -0.0990891178849438, 0.0990891178849393, 
                                 7.105427357601e-15, -1.75415237890775e-14, -7.54951656745106e-15, 
                                 -1.77635683940025e-15, 4.88498130835069e-15, 3.5527136788005e-15, 
                                 -0.099089117884942, 0.0990891178849376, -4.44089209850063e-15, 
                                 2.48689957516035e-14, -0.0990891178849393, 0.0990891178849416, 
                                 -1.77635683940025e-15, -7.21644966006352e-15, -1.02140518265514e-14, 
                                 -1.77635683940025e-15, 3.99680288865056e-15, 0, -0.0990891178849447, 
                                 0.0990891178849447, 0, 1.4432899320127e-15), Dim = c(30L, 1L), 
              Dimnames = list(NULL, NULL), factors = list())

I specify my objective function:

myobj <- function(x){
  return(-x)
}

I specify my constraint function. I have a total of 33 constraints. The last 3 constraints specified in constr2 require the solution parameters to be greater than or equal to the values in myvec.

myconstr <- function(x){
  constr1 <- as.numeric(as.matrix(cbind(mymat, -myvec2)) %*% c(x, 1))
  # Must be bigger than myvec
  constr2 <- x - myvec
  return(c(constr1, constr2))
}

Now I run the nsga2 function:

library(mco)
out <- nsga2(myobj, idim = 3, odim = 3,
                   lower.bounds = c(12.36988, -25.19478, -42.12804), 
                   upper.bounds = c(11.30719, 25.79210, 40.74460),
                   constraints = myconstr, cdim = 33)
> out$par[1, ] >= myvec
[1]  TRUE  TRUE FALSE

However, when I double check whether the solutions in out are larger than myvec, it turns out that the solutions are not satisfying the constraint that they have to be larger than myvec. Why is this the case?

Answer 1

To force the solution parameters to be greater than or equal to the values in myvec, use the lower.bounds of nsga2 instead.

You can combine your original lower.bound with myvec by getting the pairwise maximum:

pmax(myvec,c(-12.36988, -25.19478, -42.12804))

So, you get:

library(mco)

myconstr <- function(x){
  return(as.numeric(as.matrix(cbind(mymat, -myvec2)) %*% c(x, 1)))
}

out <- nsga2(myobj, idim = 3, odim = 3,
                   lower.bounds = pmax(myvec,c(-12.36988, -25.19478, -42.12804)), 
                   upper.bounds = c(11.30719, 25.79210, 40.74460),
                   constraints = myconstr, cdim = 30)
out$par[1, ] >= myvec
# [1] TRUE TRUE TRUE