2013-05-24 2 views
10

와면으로 기재된 영역 :R - 플롯 I는 다음의 부등식으로 설명 다면체를 플롯 할 RGL

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을 시도한 것이므로 코드가 다소 엉망입니다. 개선하고 싶다면 자유롭게하십시오.)

Planes

# 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 개의 변수에 대해서만 정의되지만 3 개가 있습니다. 여기서 어떻게 전진 할 수 있습니까?

+0

당신은'outer'에서 어떤 결과를 기대합니까? 'surface3d'에 대한 호출에서 어떻게 그것을 사용하려고합니까? 제발 좀 더 자세히 설명해주세요. – krlmlr

+0

예, 저는 outer를 사용하여 함수 표면을 만들고 싶습니다. – cjena

답변

5

당신은 (때문에 다른 불평등의, 교차의 일부는 다면체 외부 : 당신은뿐만 아니라 사람들을 확인해야) 한 번 에 비행기 (3)을 교차하여 다면체의 정점을 계산할 수 있습니다.

일단 정점이 있으면 연결을 시도 할 수 있습니다. 경계에있는 것을 식별하려면 세그먼트 중간 인 을 가져 와서 어떤 부등식이 평등으로 만족하는지 확인하십시오.

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

polyhedron wireframe

+0

이것은 정말 아름다운 해결책입니다! 고맙습니다. – cjena