在Mathematica中绘制线性不等式

Question

我有3个变量的线性不等式系统，我想绘制这些区域。理想情况下，我希望看起来像PolyhedronData中的对象。我试过RegionPlot3D，但结果是视觉差，过于多边形重实时旋转

这里就是我的意思是，下面的代码生成10套线性约束，并绘制他们

 
randomCons := Module[{}, 
    hadamard = KroneckerProduct @@ Table[{{1, 1}, {1, -1}}, {3}]; 
    invHad = Inverse[hadamard]; 
    vs = Range[8]; 
    m = mm /@ vs; 
    sectionAnchors = Subsets[vs, {1, 7}]; 
    randomSection := 
    Mean[hadamard[[#]] & /@ #] & /@ 
    Prepend[RandomChoice[sectionAnchors, 3], vs]; {p0, p1, p2, p3} = 
    randomSection; 
    section = 
    Thread[m -> 
     p0 + {x, y, z}.Orthogonalize[{p1 - p0, p2 - p0, p3 - p0}]]; 
    And @@ Thread[invHad.m >= 0 /. section] 
    ]; 
Table[RegionPlot3D @@ {randomCons, {x, -3, 3}, {y, -3, 3}, {z, -3, 
    3}}, {10}] 

任何建议？

更新：结合下面的建议，这里的版本，我结束了使用绘制线性不等式

 
(* Plots feasible region of a linear program in 3 variables, \ 
specified as cons[[1]]>=0,cons[[2]]>=0,... 
Each element of cons must \ 
be an expression of variables x,y,z only *) 

plotFeasible3D[cons_] := 
Module[{maxVerts = 20, vcons, vertCons, polyCons}, 
    (* find intersections of all triples of planes and get rid of \ 
intersections that aren't points *) 

    vcons = Thread[# == 0] & /@ Subsets[cons, {3}]; 
    vcons = Select[vcons, Length[Reduce[#]] == 3 &]; 
    (* Combine vertex constraints with inequality constraints and find \ 
up to maxVerts feasible points *) 
    vertCons = Or @@ (And @@@ vcons); 
    polyCons = And @@ Thread[cons >= 0]; 
    verts = {x, y, z} /. 
    FindInstance[polyCons && vertCons, {x, y, z}, maxVerts]; 
    ComputationalGeometry`Methods`ConvexHull3D[verts, 
    Graphics`Mesh`FlatFaces -> False] 
    ] 

代码的系统的可行域测试

 
randomCons := Module[{}, 
    hadamard = KroneckerProduct @@ Table[{{1, 1}, {1, -1}}, {3}]; 
    invHad = Inverse[hadamard]; 
    vs = Range[8]; 
    m = mm /@ vs; 
    sectionAnchors = Subsets[vs, {1, 7}]; 
    randomSection := 
    Mean[hadamard[[#]] & /@ #] & /@ 
    Prepend[RandomChoice[sectionAnchors, 3], vs]; {p0, p1, p2, p3} = 
    randomSection; 
    section = 
    Thread[m -> 
     p0 + {x, y, z}.Orthogonalize[{p1 - p0, p2 - p0, p3 - p0}]]; 
    And @@ Thread[invHad.m >= 0 /. section] 
    ]; 
Table[plotFeasible3D[List @@ randomCons[[All, 1]]], {50}]; 

Answer 1

这里是一个小程序，似乎做你想要什么：

rstatic = randomCons;     (* Call your function *) 
randeq = rstatic /. x_ >= y_ -> x == y; (* make a set of plane equations 
              replacing the inequalities by == *) 

eqset = Subsets[randeq, {3}];   (* Make all possible subsets of 3 planes *) 

             (* Now find the vertex candidates 
              Solving the sets of three equations *) 
vertexcandidates =  
    Flatten[Table[Solve[eqset[[i]]], {i, Length[eqset]}], 1]; 

             (* Now select those candidates 
              satisfying all the original equations *)   
vertex = Union[Select[vertexcandidates, rstatic /. # &]]; 

             (* Now use an UNDOCUMENTED Mathematica 
              function to plot the surface *) 

gr1 = ComputationalGeometry`Methods`ConvexHull3D[{x, y, z} /. vertex]; 
             (* Your plot follows *) 
gr2 = RegionPlot3D[rstatic, 
     {x, -3, 3}, {y, -3, 3}, {z, -3, 3}, 
     PerformanceGoal -> "Quality", PlotPoints -> 50] 

Show[gr1,gr2] (*Show both Graphs superposed *) 

结果是：

缺点：无证功能不健全。当面对的不是一个三角形，它会显示一个三角：

编辑

有一个选项，以摆脱犯规三角

ComputationalGeometry`Methods`ConvexHull3D[{x, y, z} /. vertex, 
           Graphics`Mesh`FlatFaces -> False] 

做的魔法。示例：

编辑2

正如我在评论中提到，这里有两套由您randomCons

{{x -> Sqrt[3/5]}, 
    {x -> -Sqrt[(5/3)] + Sqrt[2/3] y}, 
    {x -> -Sqrt[(5/3)], y -> 0}, 
    {y -> -Sqrt[(2/5)], x -> Sqrt[3/5]}, 
    {y -> 4 Sqrt[2/5], x -> Sqrt[3/5]} 
} 

和

{{x -> -Sqrt[(5/3)] + (2 z)/Sqrt[11]}, 
{x -> Sqrt[3/5], z -> 0}, 
{x -> -Sqrt[(5/3)], z -> 0}, 
{x -> -(13/Sqrt[15]), z -> -4 Sqrt[11/15]}, 
{x -> -(1/Sqrt[15]), z -> 2 Sqrt[11/15]}, 
{x -> 17/(3 Sqrt[15]), z -> -((4 Sqrt[11/15])/3)} 
} 

产生退化的顶点

还是想看看如何与那些...

编辑3

此代码轻轻应付不一般不够充分的问题，但同时也使气缸degenerancy问题为您的样本数据生成。我使用了这样的事实：致病病例总是与其轴线平行于其中一个坐标轴的圆柱体，然后使用RegionPlot3D绘制它们。 我不知道这是否会成为您的一般情况:(有用。

For[i = 1, i <= 160, i++, 
rstatic = randomCons; 
r[i] = rstatic; 
s1 = Reduce[r[i], {x, y, z}] /. {x -> var1, y -> var2, z -> var3}; 
s2 = Union[StringCases[ToString[FullForm[s1]], "var" ~~ DigitCharacter]]; 

If [Dimensions@s2 == {3}, 

    (randeq = rstatic /. x_ >= y_ -> x == y; 
    eqset = Subsets[randeq, {3}]; 
    vertexcandidates = Flatten[Table[Solve[eqset[[i]]], {i, Length[eqset]}], 1]; 
    vertex = Union[Select[vertexcandidates, rstatic /. # &]]; 
    a[i] = ComputationalGeometry`Methods`ConvexHull3D[{x, y, z} /. vertex, 
      Graphics`Mesh`FlatFaces -> False, Axes -> False, PlotLabel -> i]) 
    , 

    a[i] = RegionPlot3D[s1, {var1, -2, 2}, {var2, -2, 2}, {var3, -2, 2}, 
      Axes -> False, PerformanceGoal -> "Quality", PlotPoints -> 50, 
      PlotLabel -> i, PlotStyle -> Directive[Yellow, Opacity[0.5]], 
      Mesh -> None] 
    ]; 
] 

GraphicsGrid[Table[{a[i], a[i + 1], a[i + 2]}, {i, 1, 160, 4}]] 

Here你可以找到生成的输出的图像中，简例（所有气缸）是透明的黄色

HTH

Answer 2

A选自三重你的一组不等式通常将决定通过求解相应的三元组方程而获得的一个点。我相信你想要这个点的凸包。你可以像这样产生。

cons = randomCons; (* Your function *) 
eqs = Apply[Equal, List @@@ Subsets[cons, {3}], {2}]; 
sols = Flatten[{x, y, z} /. Table[Solve[eq, {x, y, z}], {eq, eqs}], 1]; 
pts = Select[sols, And @@ (NumericQ /@ #) &]; 
ComputationalGeometry`Methods`ConvexHull3D[pts] 

当然，一些三胞胎实际上可能是欠定的，导致线或埃文整个飞机。因此，代码将在这些情况下发出投诉。

这似乎在我尝试的几个随机病例中有效，但是，正如Yaro指出的那样，它根本不起作用。下面的图片将说明为什么：

{p0, p1, p2, 
    p3} = {{1, 0, 0, 0, 0, 0, 0, 0}, {1, 1/2, -(1/2), 0, -(1/2), 0, 
    0, -(1/2)}, {1, 0, 1/2, 1/2, 0, 0, -(1/2), 1/2}, {1, -(1/2), 1/2, 
    0, -(1/2), 0, 0, -(1/2)}}; 
hadamard = KroneckerProduct @@ Table[{{1, 1}, {1, -1}}, {3}]; 
invHad = Inverse[hadamard]; 
vs = Range[8]; 
m = mm /@ vs; 
section = 
    Thread[m -> 
    p0 + {x, y, z}.Orthogonalize[{p1 - p0, p2 - p0, p3 - p0}]]; 
cons = And @@ Thread[invHad.m >= 0 /. section]; 
eqs = Apply[Equal, List @@@ Subsets[cons, {3}], {2}]; 
sols = Flatten[{x, y, z} /. Table[Solve[eq, {x, y, z}], {eq, eqs}], 
    1]; // Quiet 
pts = Select[sols, And @@ (NumericQ /@ #) &]; 
ptPic = Graphics3D[{PointSize[Large], Point[pts]}]; 
regionPic = 
    RegionPlot3D[cons, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}, 
    PlotPoints -> 40]; 
Show[{regionPic, ptPic}] 

因此，有些点最终被其他约束定义的平面切断。这里有一个（我肯定非常低效）找到你想要的东西。

regionPts = regionPic[[1, 1]]; 
nf = Nearest[regionPts]; 
trimmedPts = Select[pts, Norm[# - nf[#][[1]]] < 0.2 &]; 
trimmedPtPic = Graphics3D[{PointSize[Large], Point[trimmedPts]}]; 
Show[{regionPic, trimmedPtPic}] 

因此，你可以使用trimmedPts的凸包。这最终取决于RegionPlot的结果，您可能需要提高PlotPoints的值以使其更可靠。

使用Google搜索可以发现线性规划中可行性区域的概念。这似乎正是你所追求的。

马克

Answer 3

看到所有的以前的答案;！有什么不对使用内置的功能 RegionPlot3D，如

RegionPlot3D[ 2*y+3*z <= 5 && x+y+2*z <= 4 && x+2*y+3*z <= 7 && 
       x >= 0 && y >= 0 && z >= 0, 
      {x, 0, 4}, {y, 0, 5/2}, {z, 0, 5/3} ]