The purpose of this notebook is to see if reasonable initial solutions to a problem articulated on OR Stack Exchange can be found via a genetic algorithm.
The problem is simple to articulate: given a matrix \(W\), find a nonzero vector \(x\) that maximizes the number of nonnegative components of \(Wx\).
First, we load libraries.
library(GA) # used to build and run a genetic algorithm
library(tictoc) # used for timing
library(magrittr) # provides a pipe operator (for more readable code)
library(ompr) # used to build a MIP model
library(ompr.roi) # along with the next two libraries, used to
library(ROI) # feed the model to CPLEX
library(ROI.plugin.cplex)We create a random problem instance. The original question looks at dimensions of 700 x 100, which we scale down to 105 x 15 so that CPLEX can find an optimal solution fairly quickly.
# Set the problem dimensions.
m <- 105
n <- 15
# Seed the random number generator.
set.seed(48823)
# Generate a random matrix W.
W <- matrix(data = runif(m * n, min = -10, max = 10), nrow = m, ncol = n)We need a function to compute the number of nonnegative entries in a vector, which will be the fitness function for the GA.
nonneg <- function(x) {
sum(W %*% x >= 0)
}We will use an “island” GA model, although a standard GA may work just as well. (The island model is intended to make premature convergence less likely). The domain for solutions is arbitrarily set to \([-10, 10]^n\). Other than population size, we will stick with the default values for parameters.
# Start timing.
tic("Running the GA:")
# Set up and run the GA.
result <- gaisl(type = "real-valued", fitness = nonneg, lower = rep(-10, n), upper = rep(10, n), popSize = 400, monitor = FALSE)
# Display the execution time.
toc()Running the GA:: 13.14 sec elapsed# Check the solution and confirm its objective value.
sol <- result@solution[1,]
cat("Solution = ", sol, "\n")Solution = 2.694808 0.542187 -1.694365 3.030605 -4.518425 0.6267654 -4.951356 3.751545 1.599644 -1.196515 -0.6927778 -0.7009321 0.4373873 4.100809 1.527456 cat("# of nonnegatives = ", nonneg(sol), "\n")# of nonnegatives = 88 To find the optimal number of nonnegative components, we run the MIP model posed in the original question (using CPLEX), with one modification. The normalization constraint is changed to use an inner product with a random vector, rather than the 1-norm of the solution.
# Set the upper limit for the "big M" constraints.
U <- n * max(abs(W))
# Generate a normalized random vector for the normalization constraint.
r <- runif(n)
r <- r / norm(r, "2")
# Create the MIP model.
mip <- MIPModel() %>%
# Add the x variables.
add_variable(x.var[i], i = 1:n, type = "continuous") %>%
# Add the binary a variables.
add_variable(a.var[i], i = 1:m, type = "binary") %>%
# Add auxiliary variables for the product W * x.
add_variable(wx[i], i = 1:m, type = "continuous") %>%
# Maximize the sum of the a variables.
set_objective(sum_expr(a.var[i], i = 1:m), sense = "max") %>%
# Define wx as W * x.
add_constraint(sum_expr(W[i, j] * x.var[j], j = 1:n) - wx[i] == 0, i = 1:m) %>%
# Add the "big M" constraints tying a to nonnegativity of wx.
add_constraint(wx[i] - U * a.var[i] <= 0, i = 1:m) %>%
add_constraint(wx[i] - U * a.var[i] >= -U, i = 1:m) %>%
# Add the normalization constraint r'x = 1.
add_constraint(sum_expr(r[i] * x.var[i], i = 1:n) == 1)
adding constraints [=>----------------------------------------] 4% eta 6s
adding constraints [=>----------------------------------------] 5% eta 5s
adding constraints [=>----------------------------------------] 6% eta 4s
adding constraints [==>---------------------------------------] 7% eta 4s
adding constraints [==>---------------------------------------] 8% eta 3s
adding constraints [===>--------------------------------------] 9% eta 3s
adding constraints [===>--------------------------------------] 10% eta 3s
adding constraints [====>-------------------------------------] 11% eta 3s
adding constraints [====>-------------------------------------] 12% eta 2s
adding constraints [=====>------------------------------------] 13% eta 2s
adding constraints [=====>------------------------------------] 14% eta 2s
adding constraints [=====>------------------------------------] 15% eta 2s
adding constraints [======>-----------------------------------] 16% eta 2s
adding constraints [======>-----------------------------------] 17% eta 2s
adding constraints [=======>----------------------------------] 18% eta 2s
adding constraints [=======>----------------------------------] 19% eta 2s
adding constraints [=======>----------------------------------] 20% eta 2s
adding constraints [========>---------------------------------] 21% eta 2s
adding constraints [========>---------------------------------] 22% eta 2s
adding constraints [=========>--------------------------------] 23% eta 2s
adding constraints [=========>--------------------------------] 24% eta 2s
adding constraints [=========>--------------------------------] 25% eta 2s
adding constraints [==========>-------------------------------] 26% eta 2s
adding constraints [==========>-------------------------------] 27% eta 2s
adding constraints [===========>------------------------------] 28% eta 2s
adding constraints [===========>------------------------------] 29% eta 2s
adding constraints [===========>------------------------------] 30% eta 1s
adding constraints [============>-----------------------------] 30% eta 1s
adding constraints [============>-----------------------------] 31% eta 1s
adding constraints [=============>----------------------------] 32% eta 1s
adding constraints [=============>----------------------------] 33% eta 1s
adding constraints [=============>----------------------------] 34% eta 1s
adding constraints [==============>---------------------------] 35% eta 1s
adding constraints [==============>---------------------------] 36% eta 1s
adding constraints [===============>--------------------------] 37% eta 1s
adding constraints [===============>--------------------------] 38% eta 1s
adding constraints [===============>--------------------------] 39% eta 1s
adding constraints [================>-------------------------] 40% eta 1s
adding constraints [================>-------------------------] 41% eta 1s
adding constraints [=================>------------------------] 42% eta 1s
adding constraints [=================>------------------------] 43% eta 1s
adding constraints [=================>------------------------] 44% eta 1s
adding constraints [==================>-----------------------] 45% eta 1s
adding constraints [==================>-----------------------] 46% eta 1s
adding constraints [===================>----------------------] 47% eta 1s
adding constraints [===================>----------------------] 48% eta 1s
adding constraints [===================>----------------------] 49% eta 1s
adding constraints [====================>---------------------] 50% eta 1s
adding constraints [=====================>--------------------] 51% eta 1s
adding constraints [=====================>--------------------] 52% eta 1s
adding constraints [=====================>--------------------] 53% eta 1s
adding constraints [======================>-------------------] 54% eta 1s
adding constraints [======================>-------------------] 55% eta 1s
adding constraints [=======================>------------------] 56% eta 1s
adding constraints [=======================>------------------] 57% eta 1s
adding constraints [=======================>------------------] 58% eta 1s
adding constraints [========================>-----------------] 59% eta 1s
adding constraints [========================>-----------------] 60% eta 1s
adding constraints [=========================>----------------] 61% eta 1s
adding constraints [=========================>----------------] 62% eta 1s
adding constraints [=========================>----------------] 63% eta 1s
adding constraints [==========================>---------------] 64% eta 1s
adding constraints [==========================>---------------] 65% eta 1s
adding constraints [===========================>--------------] 66% eta 1s
adding constraints [===========================>--------------] 67% eta 1s
adding constraints [===========================>--------------] 68% eta 1s
adding constraints [============================>-------------] 69% eta 1s
adding constraints [============================>-------------] 70% eta 1s
adding constraints [=============================>------------] 70% eta 1s
adding constraints [=============================>------------] 71% eta 0s
adding constraints [=============================>------------] 72% eta 0s
adding constraints [==============================>-----------] 73% eta 0s
adding constraints [==============================>-----------] 74% eta 0s
adding constraints [===============================>----------] 75% eta 0s
adding constraints [===============================>----------] 76% eta 0s
adding constraints [===============================>----------] 77% eta 0s
adding constraints [================================>---------] 78% eta 0s
adding constraints [================================>---------] 79% eta 0s
adding constraints [=================================>--------] 80% eta 0s
adding constraints [=================================>--------] 81% eta 0s
adding constraints [=================================>--------] 82% eta 0s
adding constraints [==================================>-------] 83% eta 0s
adding constraints [==================================>-------] 84% eta 0s
adding constraints [===================================>------] 85% eta 0s
adding constraints [===================================>------] 86% eta 0s
adding constraints [===================================>------] 87% eta 0s
adding constraints [====================================>-----] 88% eta 0s
adding constraints [====================================>-----] 89% eta 0s
adding constraints [=====================================>----] 90% eta 0s
adding constraints [=====================================>----] 91% eta 0s
adding constraints [======================================>---] 92% eta 0s
adding constraints [======================================>---] 93% eta 0s
adding constraints [=======================================>--] 94% eta 0s
adding constraints [=======================================>--] 95% eta 0s
adding constraints [=======================================>--] 96% eta 0s
adding constraints [========================================>-] 97% eta 0s
adding constraints [========================================>-] 98% eta 0s
adding constraints [=========================================>] 99% eta 0s
adding constraints [==========================================] 100% eta 0s
# Start timing the solver.
tic("Solving MIP")
# Solve the model.
result <- mip %>% solve_model(with_ROI(solver = "cplex"))CPLEX environment opened
Closed CPLEX environment# Report the run time.
toc()Solving MIP: 35.866 sec elapsed# Report the final solver status and objective value.
cat("Solver status = ", result$status, "\n")Solver status = optimal cat("Optimal objective value = ", result$objective_value, "\n")Optimal objective value = 93