Построение линейных неравенств в системе Mathematica

Вот небольшая программа, которая, кажется, делает то, что вы хотите:

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

Этот код не является достаточно общим для полной задачи, но устраняет проблему вырождения цилиндра для вашего генератора демонстрационных данных. Я использовал тот факт, что патогенные случаи всегда представляют собой цилиндры с осью, параллельной одной из осей координат, а затем использовал 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}]]

Здесь вы можете найти изображение сгенерированного вывода, выродившиеся корпуса (все цилиндры) прозрачно-желтого цвета

ХТХ!

Тройка, выбранная из вашего набора неравенств, обычно определяет точку, полученную путем решения соответствующей тройки уравнений. Я считаю, что вам нужна выпуклая оболочка этого набора точек. Вы можете сгенерировать это так.

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]

Конечно, некоторые триплеты на самом деле могут быть недоопределены и вести к линиям или даже к целой плоскости. Таким образом, в таких случаях кодекс выдает жалобу.

Это работало в нескольких случайных случаях, которые я пробовал, но, как указывает Яро, это не работает во всех случаях. Следующая картинка проиллюстрирует, почему именно:

{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 немного раскрывает концепцию области выполнимости в линейном программировании. Кажется, это именно то, что вам нужно.

отметка

Видя все предыдущие ответы; что не так с использованием встроенной функции 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} ]