使用rgl绘制由平面描述的区域

Question

我想绘制一个由以下不等式描述的多面体：

3*x+5*y+9*z<=500
4*x+5*z<=350
2*y+3*z<=150

x,y,z>=0

这是一个线性规划问题。其目标函数为：

4*x+3*y+6*z

这个程序的可行区域是一个多面体。 我可以将不等式绘制为平面，这应该可以描述多面体。 （请注意，这是我第一次尝试使用rgl，所以代码有点乱。如果您想改进它，请随意这样做）：

# setup
x <- seq(0,9,length=20)*seq(0,9,length=20)
y <- x
t <- x


f1 <- function(x,y){y=70-0.8*x}
z1 <- outer(x,y,f1)

f2 <- function(x,y){500/9-x/3-(5*y)/9}
z2 <- outer(x,y,f2)

f3 <- function(x,y){t=50-(2*y)/3}
z3 <- outer(x,y,f3)

# plot planes with rgl
uM = matrix(c(0.72428817, 0.03278469, -0.68134511, 0,
              -0.6786808, 0.0555667, -0.7267077, 0,
              0.01567543, 0.99948466, 0.05903265, 0,
              0, 0, 0, 1),
            4, 4)
library(rgl)
open3d(userMatrix = uM, windowRect = c(0, 0, 400, 400))
rgl.pop("lights")
light3d(diffuse='white',theta=0,phi=20) 
light3d(diffuse="gray10", specular="gray25")
rgl.light(theta = 0, phi = 0, viewpoint.rel = TRUE, ambient = "#FFFFFF", 
           diffuse = "#FFFFFF", specular = "#FFFFFF", x=30, y=30, z=40)
rgl.light(theta = 0, phi = 0, viewpoint.rel = TRUE, ambient = "#FFFFFF", 
          diffuse = "#FFFFFF", specular = "#FFFFFF", x=0, y=0, z=0)
bg3d("white")
material3d(col="white")
persp3d(x,y,z3,  
        xlim=c(0,100), ylim=c(0,100), zlim=c(0,100),
        xlab='x', ylab='y', zlab='z', 
        col='lightblue',
        ltheta=100, shade=0, ticktype = "simple")
surface3d(x, y, z2, col='orange', alpha=1)
surface3d(t, y, z1, col='pink', alpha=1, smooth=TRUE)

现在我想绘制由这些平面描述的区域。

x,y,z>=0.

但我不知道该怎么做。我试着这样做：

x <- seq(0,9,length=20)*seq(0,9,length=20)
y <- x
z <- x

f4 <- function(x,y,t){
  cond1 <- 3*x+5*y+9*z<=500
  cond2 <- 4*x+5*z<=350
  cond3 <- 2*y+3*z<=150

  ifelse(cond1, 3*x+5*y+9*z,
         ifelse(cond2, 4*x+5*z,
                ifelse(cond3, 2*y+3*z,0)))
}

f4(x,y,z)
z4 <- outer(x,y,z,f4) # ERROR

但是这就是我卡住的地方。outer()仅为2个变量定义，但我有三个。我该如何继续？

- cjena

---

你期望从 outer 得到什么输出？你打算如何使用它——在调用 surface3d 时？请详细说明。 - krlmlr

是的，我想使用outer函数来创建一个函数的表面。 - cjena

Answer 1

你可以通过每次相交三个平面来计算多面体的顶点 (有些交点在多面体外，因为其他不等式的原因： 你也需要检查这些)。

一旦你有了顶点，你可以尝试连接它们。 为了确定哪些在边界上，你可以取线段的中点， 并检查是否满足任何不等式作为等式。

# Write the inequalities as: planes %*% c(x,y,z,1) <= 0
planes <- matrix( c(
  3, 5, 9, -500,
  4, 0, 5, -350,
  0, 2, 3, -150,
  -1, 0, 0, 0,
  0, -1, 0, 0,
  0, 0, -1, 0
), nc = 4, byrow = TRUE )

# Compute the vertices
n <- nrow(planes)
vertices <- NULL
for( i in 1:n )
  for( j in 1:n)
    for( k in 1:n )
      if( i < j && j < k ) try( { 
        # Intersection of the planes i, j, k
        vertex <- solve(planes[c(i,j,k),-4], -planes[c(i,j,k),4] )
        # Check that it is indeed in the polyhedron
        if( all( planes %*% c(vertex,1) <= 1e-6 ) ) {
          print(vertex)
          vertices <- rbind( vertices, vertex )
        }
      } )

# For each pair of points, check if the segment is on the boundary, and draw it
library(rgl)
open3d()
m <- nrow(vertices)
for( i in 1:m )
  for( j in 1:m )
    if( i < j ) {
      # Middle of the segment
      p <- .5 * vertices[i,] + .5 * vertices[j,]
      # Check if it is at the intersection of two planes
      if( sum( abs( planes %*% c(p,1) ) < 1e-6 ) >= 2 )
        segments3d(vertices[c(i,j),])
    }

- Vincent Zoonekynd

---

这是一个非常优美的解决方案！非常感谢你。 - cjena