Solutions to a system of inequalities in R

Question

Suppose I have the following system of inequalities:

  -2x + y <= -3
1.25x + y <= 2.5
        y >= -3

I want to find multiple tuples of (x, y) that satisfy the above inequalities.

library(Rglpk)

obj <- numeric(2)
mat <- matrix(c(-2, 1, 1.25, 1, 0, 1), nrow = 3)
dir <- c("<=", "<=", ">=")
rhs <- c(-3, 2.5, -3)

Rglpk_solve_LP(obj = obj, mat = mat, dir = dir, rhs = rhs)

Using the above code only seems to return 1 possible solution tuple (1.5, 0). Is possible to return other solution tuples?

Edit: Based on the comments, I would be interested to learn if there are any functions that could help me find the corner points.

Do you have a typo in the third element of dir? Shouldn’t it be >= ? — Limey, May 29 '21 at 05:44
As there are likely infinitely many solutions, it is difficult to enumerate them all. When restricting to integer solutions or corner points things are at least theoretically simpler. Some solvers have facilities for generating some or all different integer solutions. — Erwin Kalvelagen, May 29 '21 at 12:58
@ErwinKalvelagen Makes sense. Would you happen to have recommendations on which solvers could do that? — Adrian, May 29 '21 at 14:30
High-end MIP solvers have something called a "solution pool". Constraint programming solvers can also generate multiple solutions. — Erwin Kalvelagen, May 29 '21 at 15:12
Some solvers have R interfaces, others you can use via Python (e.g. using reticulate). — Erwin Kalvelagen, May 30 '21 at 17:09

score 3 · Accepted Answer · answered May 31 '21 at 20:03

Actually to understand the possible answers for the given question we can try to solve the system of inequalities graphically.

There was a nice answer concerning plotting of inequations in R at stackowerflow. Using the given aproach we can plot the following graph:

library(ggplot2)

fun1 <- function(x) 2*x - 3        # this is the same as -2x + y <= -3
fun2 <- function(x) -1.25*x + 2.5  # 1.25x + y <= 2.5
fun3 <- function(x) -3             # y >= -3
x1 = seq(-1,5, by = 1/16)
mydf = data.frame(x1, y1=fun1(x1), y2=fun2(x1),y3= fun3(x1))
mydf <-  transform(mydf, z = pmax(y3,pmin(y1,y2)))
ggplot(mydf, aes(x = x1)) + 
  geom_line(aes(y = y1), colour = 'blue') +
  geom_line(aes(y = y2), colour = 'green') +
  geom_line(aes(y = y3), colour = 'red') +
  geom_ribbon(aes(ymin=y3,ymax = z), fill = 'gray60')

All the possible (infinite by number) tuples lie inside the gray triangle. The vertexes can be found using the following code.

obj <- numeric(2)
mat <- matrix(c(-2, 1.25, 1, 1), nrow = 2)
rhs <- matrix(c(-3, 2.5), nrow = 2)

aPoint <- solve(mat, rhs)

mat <- matrix(c(-2, 0, 1, 1), nrow = 2)
rhs <- matrix(c(-3, -3), nrow = 2)

bPoint <- solve(mat, rhs)

mat <- matrix(c(1.25, 0, 1, 1), nrow = 2)
rhs <- matrix(c(2.5, -3), nrow = 2)

cPoint <- solve(mat, rhs)

Note the order of arguments of matrices.

And you get the coordinates:

> aPoint
          [,1]
[1,] 1.6923077
[2,] 0.3846154
> bPoint
     [,1]
[1,]    0
[2,]   -3
> cPoint
     [,1]
[1,]  4.4
[2,] -3.0

ThomasIsCoding · Answer 2 · 2021-06-03T15:53:37.073

All the codes below are with base R only (no need library(Rglpk))

1. Corner Points

If you want to get all the corner points, here is one option

A <- matrix(c(-2, 1.25, 0, 1, 1, -1), nrow = 3)

b <- c(-3, 2.5, 3)

# we use `det` to check if the coefficient matrix is singular. If so, we return `Inf`.
xh <-
  combn(nrow(A), 2, function(k) {
    if (det(A[k, ]) == 0) {
      rep(NA, length(k))
    } else {
      solve(A[k, ], b[k])
    }
  })

# We filter out the points that satisfy the constraint
corner_points <- t(xh[, colSums(A %*% xh <= b, na.rm = TRUE) == length(b)])

such that

> corner_points
         [,1]       [,2]
[1,] 1.692308  0.3846154
[2,] 0.000000 -3.0000000
[3,] 4.400000 -3.0000000

2. Possible Tuples

If you want to have multiple tuples, e.g., n=10, we can use Monte Carlo simulation (based on the obtained corner_points in the previous step) to select the tuples under the constraints:

xrange <- range(corner_points[, 1])
yrange <- range(corner_points[, 2])
n <- 10
res <- list()
while (length(res) < n) {
  px <- runif(1, xrange[1], xrange[2])
  py <- runif(1, yrange[1], yrange[2])
  if (all(A %*% c(px, py) <= b)) {
    res[length(res) + 1] <- list(c(px, py))
  }
}

and you will see n possible tuples in a list like below

> res
[[1]]
[1]  3.643167 -2.425809

[[2]]
[1]  2.039007 -2.174171

[[3]]
[1]  0.4990635 -2.3363637

[[4]]
[1]  0.6168402 -2.6736421

[[5]]
[1]  3.687389 -2.661733

[[6]]
[1]  3.852258 -2.704395

[[7]]
[1] 1.7571062 0.1067597

[[8]]
[1]  3.668024 -2.771307

[[9]]
[1]  2.108187 -1.365349

[[10]]
[1]  2.106528 -2.134310

Waldi · Answer 3 · 2021-05-31T21:43:27.600

First of all, the matrix representing the three equations needs a small correction, because R fills matrices column by column :

  -2x + y <= -3
1.25x + y <= 2.5
        y >= -3

mat <- matrix(c(-2, 1.25, 0, 1, 1, 1), nrow = 3
# and not : mat <- matrix(c(-2, 1, 1.25, 1, 0, 1), nrow = 3)

To get different tuples, you could modify the objective function :

obj <- numeric(2) results in an objective function 0 * x + 0 * y which is always equal to 0 and can't be maximized : the first valid x,y will be selected.

Optimization on x is achieved by using obj <- c(1,0), resulting in maximization / minimization of 1 * x + 0 * y.
Optimization on y is achieved by using obj <- c(0,1).

#setting the bounds is necessary, otherwise optimization occurs only for x>=0 and y>=0
bounds <- list(lower = list(ind = c(1L, 2L), val = c(-Inf, -Inf)),
               upper = list(ind = c(1L, 2L), val = c(Inf, Inf)))

# finding maximum x: obj = c(1,0), max = T
Rglpk_solve_LP(obj = c(10,0), mat = mat, dir = dir, rhs = rhs,bound=bounds, max = T)$solution
# [1] 4.4 -3.0

# finding minimum x: obj = c(1,0), max = F
Rglpk_solve_LP(obj = c(10,0), mat = mat, dir = dir, rhs = rhs,bound=bounds, max = F)$solution
#[1]  0 -3

# finding maximum y: obj = c(0,1), max = T
Rglpk_solve_LP(obj = c(0,1), mat = mat, dir = dir, rhs = rhs,bound=bounds, max = T)$solution
#[1] 1.6923077 0.3846154

Solutions to a system of inequalities in R

3 Answers3

1. Corner Points

2. Possible Tuples