R - Plot eine Region von Flugzeugen mit rgl beschrieben

Question

ich einen Polyeder dargestellt werden soll, die durch die folgenden Ungleichungen beschrieben:

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

x,y,z>=0 

Es ist ein lineares Programm. Die Zielfunktion lautet:

Das Polyeder ist die mögliche Region für dieses Programm. Ich bin in der Lage, die Ungleichungen als Ebenen zu plotten, die das Polyeder beschreiben sollten (Beachten Sie, dass dies mein erster Versuch mit rgl ist, so ist der Code irgendwie chaotisch. Wenn Sie es verbessern möchten, zögern Sie bitte nicht):

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

Jetzt möchte ich die Region zeichnen, die von den Ebenen mit

x,y,z>=0. 

Aber ich weiß nicht beschrieben, wie es zu tun. Ich habe versucht, es so zu machen:

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 

Aber das ist der Punkt, wo ich feststecke. outer() ist nur für 2 Variablen definiert, aber ich habe drei. Wie kann ich von hier weiterziehen?

Answer 1

Sie können die Ecken des Polyeders berechnen, indem die Ebenen 3 zu einer Zeit, sich schneid (einige der Kreuzungen außerhalb der Polyeder sind, weil andere Ungleichheiten: Sie haben diese ebenfalls zu überprüfen).

Sobald Sie die Scheitelpunkte haben, können Sie versuchen, sie zu verbinden. Um zu identifizieren, welche an der Grenze sind, können Sie die Mitte des Segments nehmen und prüfen, ob eine Ungleichheit als eine Gleichheit erfüllt ist.

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