R - Plot a region described by planes with rgl

Question

I want to plot a polyhedron, which is described by the following inequalities:

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

x,y,z>=0

It is a linear program. The objective function is:

4*x+3*y+6*z

The polyhedron is the feasible region for this program. I am able to plot the inequalities as planes, which should describe the polyhedron (Note that this is my first try with rgl, so the code is kinda messy. if you want to improve it, please feel free to do so):

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

Now I want to plot the region that is described by the planes with

x,y,z>=0.

But I don't know how to do it. I tried to do it like this:

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

But this is the point where I'm stuck. outer() is defined only for 2 variables, but I have three. How can I move on from here?

Answer 1

You can compute the vertices of the polyhedron by intersecting the planes 3 at a time (some of the intersections are outside the polyhedron, because of other inequalities: you have to check those as well).

Once you have the vertices, you can try to connect them. To identify which are on the boundary, you can take the middle of the segment, and check if any inequality is satisfied as an equality.

# 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),])
    }