Setting constraints of a multi-objective function in nsga2 R

Question

I want to use the function nsga2 in library mco to solve a multi-objective problem and find the Pareto Frontier but I couldn't set up the constraints properly.

The objective functions are below. The context of the problem is project selection, i.e., I have five projects represented by x[1], x[2], ... x[5] and only some can be selected. For example, if project No. 1 is selected, then x[1]=1 if not selected x[1]=0 and this is true for all projects (value of x[n] is discrete, either 1 or 0). Another constraint I have is the total budget of the selected projects should be less than 100. After running the nsga2 function, the parameters in the Solution doesn't seem to be right as the parameters are not 1 or 0. Do I have the constraints wrong? How can I find optimum values of x[1] to x[5]? Thanks!

# objective functions to minimize
 ObjFun <- function (x){
   f1 <- -0.02*x[1] + 0.01*x[2] + 0.02*x[3] + -0.01*x[4] + 0.02*x[5]  
   f2 <- 0.17*x[1] + -0.08*x[2] + 0.10*x[3] + 0.09*x[4] + 0.07*x[5] 
   c(f1, f2)      }

 # The constraints 
 Constr <- function(x){ 
   100 >= 20*x[1] + 30*x[2] + 20*x[3] + 33*x[4] + 60*x[5] # Total budget >= total project costs
   x[1:5] == 1
   x[1:5] == 0      }

 library(mco)
 Solution <- nsga2(ObjFun, 5, 2, lower.bounds=c(0,0,0,0,0), upper.bounds=c(1,1,1,1,1), constraints = Constr)
 # plot(Solution)
 Solution$par

Answer 1

Since x[i] can only be 1 or 0, you are dealing with a combinatorial optimization problem, where the space you have to optimize on is discrete:

https://en.wikipedia.org/wiki/Combinatorial_optimization

In general, numerical optimization routines are constructed to work on continuous spaces (subsets of R^n). However, in your case the discrete space is small and the problem lends itself to a simple brute force approach, where you evaluate ObjFunc on all the 32 possible points. The Pareto frontier is here also discrete.

## objective functions to minimize
ObjFun <- function (x){
    f1 <- -0.02*x[1] + 0.01*x[2] + 0.02*x[3] + -0.01*x[4] + 0.02*x[5]  
    f2 <- 0.17*x[1] + -0.08*x[2] + 0.10*x[3] + 0.09*x[4] + 0.07*x[5] 
    c(f1=f1, f2=f2)
}

## space of all 32 feasible solutions 
space <- expand.grid(data.frame(matrix(0:1, nrow=2, ncol=5)))
## brute force evaluation of ObjFun on all the 32 feasible solutions
val <- sapply(data.frame(t(space)),  ObjFun)
tmp <- sol <- cbind(space, t(val))

## returns indices of all rows which are Pareto dominated
## by the i-th row
which.are.dominated <- function(i, tmp){
    s1 <- tmp$f1[i]
    s2 <- tmp$f2[i]
    with(tmp,
         which( (s1 <= f1) &
                (s2 <= f2) &
                ( (s1 < f1) |
                  (s2 < f2) )
               ))
}
## For each feasible solution i, remove all feasible solutions which are Pareto dominated by feasible solutions i
i <- 1
repeat{
    remove <- which.are.dominated(i, tmp)
    if(length(remove)>0) tmp <- tmp[-remove, ]
    if(i>=nrow(tmp)) break
    i <- i+1
}
with(sol, plot(f1, f2))
points(tmp$f1, tmp$f2, pch=20, col=2)
legend("topright", col=2, pch=20, "Pareto frontier")

References:

https://en.wikipedia.org/wiki/Multi-objective_optimization

https://en.wikipedia.org/wiki/Pareto_efficiency

P.S. I got to use a repeat statement for possibly the first time ever since I started to use R years ago...

Edit: A non-brute-force approach is to use nsga2 :D As I set it up, solutions are searched for x varying in the n-dimensional cube [0,1]^n where n is the number of projects; the algorithm produces a number of solutions (200 in my example), which you can then "discretize" to 0 or 1 using round. For larger numbers of projects, to obtain a more precise approximation to the Pareto frontier, you have to use more generations (e.g. 600). In the final plots, only a sample of costs is plotted if more than 12 projects are considered.

##n.projects <- 12
n.projects <- 50
if(n.projects>25) generations=600


set.seed(1)
vecf1 <- rnorm(n.projects)
vecf2 <- rnorm(n.projects)
vcost <- rnorm(n.projects)
n.solutions <- 200

library(mco)

ObjFun <- function (x){
    f1 <- sum(vecf1*x)
    f2 <- sum(vecf2*x)
    c(f1=f1, f2=f2)
}
Constr <- function(x){ 
    c(100 - sum(vcost*x)) # Total budget >= total project costs
}

Solution <- nsga2(ObjFun, n.projects, 2,
                  lower.bounds=rep(0,n.projects), upper.bounds=rep(1,n.projects),
                  popsize=n.solutions,  constraints = Constr, cdim=1,
                  generations=generations)

selected.project.combinations <- unique(round(Solution$par))
selected.project.combinations.costs <- sapply(data.frame(t(selected.project.combinations)),  ObjFun)

## final plotting of results
max.n.proj.plot <- 12
if(n.projects <= max.n.proj.plot){
    xsamp <- expand.grid(data.frame(matrix(0:1, nrow=2, ncol=n.projects)))
}else{
    xsamp <- matrix(sample(0:1, n.projects*2^max.n.proj.plot, replace=TRUE), ncol=n.projects)
}
fsamp <- sapply(data.frame(t(xsamp)), ObjFun)

par(mfrow=c(1,2))
plot(Solution)
points(fsamp[1, ], fsamp[2, ])
points(t(selected.project.combinations.costs), col=3, pch=20)
legend("bottomleft", bty="n", pch=c(20,1), col=c(3,1),
       c("Costs of optimal\nproject combinations",
         "Costs of discarded\nproject combinations"),
       y.intersp=1.8
       )

plot(t(fsamp), xlim=range(Solution$value[ ,1], fsamp[1, ]),
     ylim=range(Solution$value[ ,2], fsamp[2, ]))
points(Solution$value, col=2, pch=".")
points(t(selected.project.combinations.costs), col=3, pch=20)