R中不等式系统的解决方案

Question

假设我有以下不等式系统：

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

我想找到满足上述不等式的(x, y) 的多个元组。

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)

使用上面的代码似乎只返回 1 个可能的解决方案元组(1.5, 0)。是否可以返回其他解决方案元组？

编辑：基于 cmets，我有兴趣了解是否有任何功能可以帮助我找到角点。

Answer 1

实际上，要了解给定问题的可能答案，我们可以尝试以图形方式解决不等式系统。

在 stackowerflow 上有一个不错的 answer concerning plotting of inequations in R。使用给定的方法，我们可以绘制以下图表：

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')

所有可能的（无限数量的）元组都位于灰色三角形内。 可以使用以下代码找到顶点。

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) 

注意矩阵参数的顺序。

然后你得到坐标：

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

Answer 2

以下所有代码均以 R 为基础（不需要library(Rglpk)）

---

1。角点

如果你想得到所有的角点，这里有一个选项

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)])

这样

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

---

2。可能的元组

如果你想有多个元组，例如n=10，我们可以使用蒙特卡罗模拟（基于上一步得到的corner_points）来选择约束下的元组：

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))
  }
}

你会在下面的列表中看到n可能的元组

> 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

Answer 3

首先，表示三个方程的矩阵需要稍作修正，因为 R 逐列填充矩阵：

  -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)

要获得不同的元组，您可以修改目标函数：

obj <- numeric(2) 导致目标函数 0 * x + 0 * y 始终等于 0 且无法最大化：将选择第一个有效的 x,y。

对x 的优化是通过使用obj <- c(1,0) 实现的，从而最大化/最小化1 * x + 0 * y。
对y的优化是通过使用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