Построение линейных неравенств в Mathematica - PullRequest
4 голосов
/ 28 сентября 2010

У меня есть линейные системы неравенств в 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}];

Ответы [ 3 ]

3 голосов
/ 06 октября 2010

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

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

Результат:

alt text

Недостаток: недокументированная функцияне идеально.Если грань не является треугольником, она будет отображать триангуляцию:

alt text

Редактировать

Существует возможность избавиться отгрязная триангуляция

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

делает магию.Пример:

alt text

Редактировать 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}]]

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

HTH!

3 голосов
/ 29 сентября 2010

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

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, чтобы сделать его более надежным.

Гугление немного раскрывает концепцию области осуществимости в линейном программировании.Кажется, это именно то, что вам нужно.

Mark

1 голос
/ 04 мая 2011

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