Message Boards Message Boards

Packing multiple ellipses in a non-convex curve

enter image description here

Introduction

Packing ellipses in a non-convex oval involves two problems: keeping the ellipses inside the oval and keep the ellipses disjoint. The basic approach for keeping the ellipses inside was discussed in a previous post

http://community.wolfram.com/groups/-/m/t/1453823

This function is applied to a group of discrete points, from which an interpolation function is created for use in the optimization.
Keeping the ellipses disjoint, a simpler problem, is done by using Lagrange multipliers to create constraints, see

http://www.optimization-online.org/DB_FILE/2016/03/5348.pdf

The optimization is done as as random searches using ParametricIPOPTMinimize.

 ParametricIPOPT Code

In[1]:=  Needs["IPOPTLink`"] 

In[2]:= statusrl = {0 -> "Solve_Succeeded", 1 -> "Solved_To_Acceptable_Level", 
   2 -> "Infeasible_Problem_Detected", 
   3 -> "Search_Direction_Becomes_Too_Small", 4 -> "Diverging_Iterates", 
   5 -> "User_Requested_Stop", 
   6 -> "Feasible_Point_Found", -1 -> "Maximum_Iterations_Exceeded", -2 -> 
    "Restoration_Failed", -3 -> "Error_In_Step_Computation", -4 -> 
    "Maximum_CpuTime_Exceeded", -10 -> "Not_Enough_Degrees_Of_Freedom", -11 ->
     "Invalid_Problem_Definition", -12 -> "Invalid_Option", -13 -> 
    "Invalid_Number_Detected", -100 -> "Unrecoverable_Exception", -101 -> 
    "NonIpopt_Exception_Thrown", -102 -> "Insufficient_Memory", -199 -> 
    "Internal_Error"}; 

In[3]:=  

In[4]:= iMin[obj_, cons_List, vwb_List, nrands_Integer, seed_Integer, opts___] :=
 Block[{lecons, vars, lbubs, cfuns, clbubs, v, lb, ub, params, param, pip, 
   starts, pres, listres, status, goodres, finalres},

  lecons = cons /. {
        ( a_ == b_  ) -> LessEqual[0, a - b, 0],
        ( a_ <=  b_ ) -> LessEqual[-\[Infinity], a - b, 0],
        (a_ >= b_  ) -> LessEqual[0, a - b, \[Infinity]]
     };

  {vars, lbubs} = Transpose[vwb /. 
     {v_, lb_?NumericQ, ub_?NumericQ} :> {v, {lb, ub}}];

  If[Length[cons] == 0, {cfuns, clbubs} = { {}, {} },
   {cfuns, clbubs} = 
    Transpose[lecons /. LessEqual[lb_, v_, ub_] -> {v, {lb, ub}} ]
   ];

  params = Table[param[i], {i, Length[vwb]}];
  pip = ParametricIPOPTMinimize[obj, vars, params, lbubs, cfuns, clbubs, 
    params, "RuntimeOptions" -> {"WarningMessages" -> False}, opts];
  SeedRandom[seed];
  starts = Transpose[RandomReal[#, nrands] & /@ lbubs];
  pres = pip @@@ starts;
  listres = {IPOPTMinValue[#], IPOPTArgMin[#], IPOPTReturnCode[#]} & /@ pres;
  IPOPTDataDelete /@ pres;
  status = listres[[All, 3]];
  Print[Tally[status] /. {s_, n_Integer} :> { s /. statusrl, n}];
  goodres = Select[listres, Or[#[[3]] == 0, #[[3]] == 1] &];

  finalres = If[Length[goodres] == 0,
    listres[[Ordering[listres[[All, 1]]]]][[1]],
    goodres[[ Ordering[goodres[[All, 1]]]]][[1]]
    ];

  {finalres[[1]],
   Thread[vars -> finalres[[2]]], finalres[[3]] /. statusrl}
  ]

 Packing Code

 equation of circumscribing oval is oval[{x,y}] == 0

In[5]:= oval[{x_, y_}] = ((x - 1)^2 + y^2) ((x + 1)^2 + y^2) - (21/20)^4;

code to generate equation of an ellipse with semi-major axes "a" and "b" 

In[6]:= rottransrl[{xc_, yc_, \[Theta]_}, {x_, y_}] = 
  Thread[{x, y} -> RotationMatrix[-\[Theta]].{x - xc, y - yc}];

In[7]:= Thread[{xel, yel} =  
   DiagonalMatrix[{a, b}^-1].RotationMatrix[-\[Theta]].{x - xc, y - yc}];

In[8]:= ell[a_, b_][xc_, yc_, \[Theta]_][x_, y_] =  xel^2 + yel^2 - 1;

 eliminate \[Lambda] from Lagrange multiplier equation for minimizing size of ellipse with axes "a" and "a/2" that overlaps oval

In[9]:= eld = Eliminate[
   D[oval[{x, y}] == \[Lambda] ell[a, a/2][xc, yc, \[Theta]][x, y], {{x, 
      y}}], \[Lambda]] // FullSimplify

Out[9]= 5 y (xc + x (-2 + x xc) + xc y^2) + 
   3 (-xc y + (2 x - xc) y (x^2 + y^2) - x (-1 + x^2 + y^2) yc) Cos[
     2 \[Theta]] == 
  5 x (-1 + x^2 + y^2) yc + 
   3 (x^4 - x^3 xc + x (xc - xc y^2) - y (1 + y^2) (y - yc) + 
      x^2 (-1 + y yc)) Sin[2 \[Theta]] && a != 0

 define function using previous result without a!=0

In[10]:= el[xc_, yc_, \[Theta]_] = eld[[1]]

Out[10]= 5 y (xc + x (-2 + x xc) + xc y^2) + 
  3 (-xc y + (2 x - xc) y (x^2 + y^2) - x (-1 + x^2 + y^2) yc) Cos[
    2 \[Theta]] == 
 5 x (-1 + x^2 + y^2) yc + 
  3 (x^4 - x^3 xc + x (xc - xc y^2) - y (1 + y^2) (y - yc) + 
     x^2 (-1 + y yc)) Sin[2 \[Theta]]

 define function to find larger axis of largest ellipse centered at {xc,yc} with orientation \[Theta] that finds inside oval.  
Min is needed due to multiple solution of Lagrange multiplier equations.

In[11]:= dist[{xc_?NumericQ, 
   yc_?NumericQ, \[Theta]_?NumericQ}] := -N @ Sign[oval[{xc, yc}]]*
  Min[a /. NSolve[{a >= 0,
      oval[{x, y}] == 0,
      ell[a, a/2][xc, yc, \[Theta]][x, y] == 0,
      el[xc, yc, \[Theta]]},
     {a, x, y}, Reals]]

 evaluate the function at certain discrete values

In[12]:= AbsoluteTiming[
 Quiet[disttblflat = 
    N @ Flatten[
      Table[{xc, yc, \[Theta], dist[{xc, yc, \[Theta]}]}, {xc, -2, 2, 1/
        10}, {yc, -1, 1, 1/10}, {\[Theta], 0, 2 \[Pi], \[Pi]/10}], 2]];]

Out[12]= {797.969, Null}

 construct an interpolation function from the values at discrete points

In[13]:= Dimensions[
 disttbl = disttblflat /. {a_?NumericQ, b_, c_, d_} :> {{a, b, c}, d}]

Out[13]= {18081, 2}

In[14]:= itbl = Interpolation[disttbl, InterpolationOrder -> 2];

 area of the oval for calculating packing fraction

In[15]:= ovalarea = NIntegrate[Boole[oval[{x, y}] <= 0], {x, -2, 2}, {y, -2, 2}]

Out[15]= 2.56423

 plot of the oval

In[16]:= ovalplot = 
  ContourPlot[Evaluate[oval[{x, y}] == 0], {x, -1.5, 1.5}, {y, -1.5, 1.5}, 
   ImageSize -> Small, ContourStyle -> Red];

 forall[vars, c == d, a <= b] using Lagrange multipliers

In[17]:= lagforall[LessEqual[a_, b_], Equal[c_, d_], vars_List] :=
 {Flatten @ {LessEqual[a, b], Equal[c, d], \[Lambda] >=  0, 
     D[a - b == \[Lambda] (c - d), {vars}]} /. \[Lambda] -> \[Lambda][
     Sequence @@ vars],
   Join[vars, {\[Lambda][Sequence @@ vars]}]}

 constraint for two curves to not overlap

In[18]:=  curvesnooverlap[curve1_, {xc1_, yc1_, \[Theta]1_}, 
  curve2_, {xc2_, yc2_, \[Theta]2_}, {x_, y_}] :=
 lagforall[
  0 <= curve1[xc1, yc1, \[Theta]1][x, y],
  curve2[xc2, yc2, \[Theta]2][x, y] == 0,
  {x, y}
  ]

no overlap constraint for two ellipses 

In[19]:= noeqs[i_, j_] = 
  curvesnooverlap[ell[a, a/2], {xc[i], yc[i], \[Theta][i]}, 
   ell[a, a/2], {xc[j], yc[j], \[Theta][j]}, {xno[i, j], yno[i, j]}];

 all the constraints:  all the ellipses fit into the oval and all are disjoint

In[20]:= cons[n_Integer] := 
 Join[Table[a == itbl[xc[i], yc[i], \[Theta][i]], {i, n}], 
  Flatten[Table[Sequence @@ noeqs[i, j][[1]], {i, n - 1}, {j, i + 1, n}], 1]]

 all the variables with bounds

In[21]:= vwb[n_Integer] := 
 Join[{{a, 0, 1}}, Table[{xc[i], -1.5, 1.5}, {i, n}], 
  Table[{yc[i], -.5, .5}, {i, n}], 
  Table[{\[Theta][i], \[Pi]/2, 3 \[Pi]/2}, {i, n}], 
  Flatten[Table[{xno[i, j], -1.5, 1.5}, {i, n - 1}, {j, i + 1, n}], 1], 
  Flatten[Table[{yno[i, j], -.5, .5}, {i, n - 1}, {j, i + 1, n}], 1], 
  Flatten[Table[{\[Lambda][xno[i, j], yno[i, j]], 0, 5}, {i, n - 1}, {j, 
     i + 1, n}], 1]]

 plotting function, putting packing fraction over plot

In[22]:= plt[n_Integer, sln_] := 
 Show[ovalplot, 
  ContourPlot[
   Evaluate[Table[ell[a, a/2][xc[i], yc[i], \[Theta][i]][x, y] == 0, {i, n}] /. 
     sln[[2]]], {x, -1.5, 1.5}, {y, -1.5, 1.5}],
  PlotLabel -> ToString[ ((n*\[Pi] a a/2)/ovalarea /. sln[[2]])]
  ]

 Packing Calculation

In[23]:= plts = Range[10];

In[24]:= AbsoluteTiming[plts[[1]] = plt[1, iMin[-a, cons[1], vwb[1], 100, 0]]];

During evaluation of In[24]:= {{Maximum_Iterations_Exceeded,3},{Solve_Succeeded,97}}

In[25]:= AbsoluteTiming[
 Do[plts[[i]] = plt[i, iMin[-a, cons[i], vwb[i], 50, 0]], {i, 2, 10}]]

During evaluation of In[25]:= {{Solve_Succeeded,40},{Maximum_Iterations_Exceeded,5},{Infeasible_Problem_Detected,5}}

During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,26},{Infeasible_Problem_Detected,13},{Restoration_Failed,1},{Solve_Succeeded,10}}

During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,31},{Infeasible_Problem_Detected,12},{Solve_Succeeded,7}}

During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,36},{Solve_Succeeded,9},{Infeasible_Problem_Detected,5}}

During evaluation of In[25]:= {{Solve_Succeeded,11},{Maximum_Iterations_Exceeded,34},{Infeasible_Problem_Detected,5}}

During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,35},{Solve_Succeeded,9},{Infeasible_Problem_Detected,6}}

During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,43},{Solve_Succeeded,6},{Infeasible_Problem_Detected,1}}

During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,35},{Solve_Succeeded,8},{Infeasible_Problem_Detected,7}}

During evaluation of In[25]:= {{Maximum_Iterations_Exceeded,36},{Solve_Succeeded,7},{Infeasible_Problem_Detected,7}}

Out[25]= {3989.62, Null}

In[26]:= plts
POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract