GA optimization of 3D point cloud

I have a 3D point cloud I need to optimize with respect to maximizing the volume of its convex hull. For n points, the input to the fitness function currently looks like:

arr(3 * n) = [ X1, Y1, Z1, X2, Y2, Z2, ... Xn, Yn, Zn ]

This works just fine till I add upper and lower bounds to the Z coordinates via vectors which look like:

LB(3 * n) = [ -Inf, -Inf, 3000, ... -Inf, -Inf, 3000 ]
UB(3 * n) = [  Inf,  Inf, 9000, ...  Inf,  Inf, 9000 ]

Here, the optimization performance and solution quality both drop greatly. The convex hull of the point cloud now tends to approach a cube, instead of freely expanding in X and Y. I think this because the GA only sees them as a vector of doubles, being unable to distinguish between X, Y, and Z coordinates, having no clue of their spatial significance.

I tried to implement custom data type optimization using cells, but run into other problems, the most major of which is of how to specify upper and lower bounds (which GA insists must be double vectors). The documentation for custom data type optimization offers no insight into this either.

Another wild guess would be to create 3 subpopulations of X, Y, and Z coordinates, and optimize them within themselves with no crossover, but I'm not sure if this would be feasible.

I would be grateful if someone could point me towards a viable approach for modelling this scenario, and would be glad to provide relevant snippets of code for reference as required.

EDIT (fitness function):

   function volume = fitfun(N, credit, NUM)

     INDIVS = size(N, 1);
     volume = zeros(INDIVS, 1);
     overlapFraction = 0.2;

     parfor i = 1 : INDIVS
       inferior = false;
       N2 = reshape(N(i,:), 3, NUM).';
       [ facets, volume(i) ] = convexHull(delaunayTriangulation(N2));
       numFacets = size(facets, 1);
       numVerts = size(facets, 2);

       for f = 1 : numFacets
         for v1 = 1 : numVerts
           v2 = rem(v1, numVerts) + 1;
           p = [ facets(f,v1); facets(f,v2) ];
           n = [ N2(p(1),:); N2(p(2),:) ];
           range = [ credit(p(1)); credit(p(2)) ];
           edge = norm(n(1,:) - n(2,:));

           % Edge length constraint:
           range = range_calc(n, range, edge);

           coverage = range(1) + range(2);
           overlap = coverage - edge;
           minOverlap = edge * overlapFraction;

           if overlap < minOverlap
             volume(i) = 0;
             inferior = true;
             break
           end
         end

         if inferior == true
           break
         end
       end

       if inferior == true
         continue
       end
     end

   end

[Matt J]

I have a 3D point cloud I need to optimize with respect to maximizing the volume of its convex hull.

If you have no bounds on Xi,Yi, the maximum volume of the hull should be Inf.

[slaiyer]

There are edge length constraints for every pair of neighboring nodes that are enforced quite satisfactorily within the fitness function for the current unbounded scenario. This issue shows itself when LB and UB are introduced for the Z coordinates.

[Matt J]

You should probably post the fitness function and constraints. It might just be that there are no satisfactory solutions when Z is bounded to the interval you've indicated.

[slaiyer]

I have appended the fitness function to the original post.

It seems unlikely to me that the interval in question Z = [ 3000, 9000 ] is the issue as it is almost always a subset of the Z interval of the unbounded solutions, which are veritably optimal but for the lack of Z bounds.

The nature of the problem is such that it is guaranteed to be valid in Z = [ 0, 10000 ]

[Matt J]

I think we need to see the constraints, too.

Note also that the solution will be ambiguous unless one of the vertices is fixed at a known location. Since translating the polyhedron won't change its volume, the solution could only be unique up to a translation of all the vertices. Rotations and other isometric transformations would be issues that way as well....

[slaiyer]

Translation, rotation, and other isometric transformations on the polyhedron will not be valid as the constraints are fed with data set values, which vary with Z.

range_calc() and its callees contain way too much code to be attached here. The summary would be that it returns ranges in [ 3000, 6000 ] for both vertices of each edge. I can assure you that they are very much valid in Z = [ 3000, 9000 ], and that ranges are independent of X and Y.

It is just a question of squashing the unbounded solution into Z = [ 3000, 9000 ]

[Matt J]

It seems unlikely to me that the interval in question Z = [ 3000, 9000 ] is the issue as it is almost always a subset of the Z interval of the unbounded solutions

If the solution without the imposed Z-bounds already lies in the interval [3000,9000] then

1. Why do you feel the need to impose the bounds at all?
2. What happens when you put that solution in the initial population and re-run with the bounds imposed?

[slaiyer]

1. The particular scenario for which I'm deriving the solution has these bounds on the Z domain.

2. When I feed the final population (5th return value) of the unbounded GA to the bounded GA, it behaves exactly like it does when the bounded GA is run from scratch: it tries to create a point cloud whose convex hull approaches a cube as it cannot differentiate between X, Y, and Z coordinates.

Just to be clear, I would want a disc-like expansion in the XY plane with only Z bounded to [ 3000, 9000], instead of all three axes being restricted despite the explicitly specified infinite bounds for X and Y.

[Matt J]

When I feed the final population (5th return value) of the unbounded GA to the bounded GA, it behaves exactly like it does when the bounded GA is run from scratch: it tries to create a point cloud whose convex hull approaches a cube as it cannot differentiate between X, Y, and Z coordinates.

I've lost track of what you're claiming. You claimed earlier (I thought) that unbounded GA already produces satisfactory results, i.e., the Z's returned by unbounded GA turn out to be already bounded to the interval [3000, 9000] even though you didn't impose that explicitly through the LB,UB arguments.

If that's true, then putting that solution in the initial population of bounded GA should cause GA to stop after 0 generations, because the solution is already optimal and satisfies all the bounds already including [3000,9000].

If that's not what you're claiming then please explain again why you think imposing bounds on Z shouldn't change the solution drastically.

[slaiyer]

Sorry, I must not have been clear in my description of the situation.

The unbounded GA is constrained only by the edge length constraints specified. It is free to expand in all directions roughly around whatever origin is implied in 'PopInitRange' (or in 'InitialPopulation' when I'm testing with an explicit population), the score of an individual being the volume of the convex hull of the point cloud it forms. It may end up occupying space in say Z = [ -4000, 6000 ] if I center the initial population around [ 0, 0, 0 ].

(What I meant was that [ 3000, 9000 ] is almost always a subset of the Z interval of the unbounded GA's solution, especially when I center the initial population around a point anywhere in Z = 6000. I mentioned this to lay to rest your concerns that there may be no satisfactory solutions in this interval.)

Since the current scenario to be dealt with has Z bounded to [ 3000, 9000 ], I center the initial population around [ 6000, 6000, 6000 ] (any X and Y values will do). The unbounded solution may now end up in Z = [ 0, 12000 ] (of which my target bounds of [ 3000, 9000 ] are clearly a subset).

To derive a bounded solution in [ 3000, 9000 ], I obviously center the initial population around some point in Z = 6000 ((LB + UB) / 2) for both the bounded (my target) and unbounded (currently reliable) GAs.

I need to find a way to restrict the otherwise unbounded solution to the interval [ 3000, 9000 ] along the Z axis, while letting it freely expand in X and Y, essentially creating a very thick pancake of sorts, if I may. I hope I have clarified my predicament.

[Matt J]

The unbounded solution may now end up in Z = [ 0, 12000 ] (of which my target bounds of [ 3000, 9000 ] are clearly a subset).

Well...there's little that we can conclude from that. We also know trivially that the unbounded solution will end up in Z=[-inf, inf] of which [3000,9000] is clearly a subset.

What exactly is defective in the solution that bounded GA gives you? Does it not satisfy all the constraints?

[slaiyer]

My intention in providing the example intervals of Z = [ -4000, 6000 ] (initial population centered around [ X_any, Y_any, 0 ]) and Z = [ 0, 12000 ] (initial population centered around [ X_any, Y_any, 6000 ]) was to convey an idea of the general length of the unbounded solution interval along the Z axis (restricted to such sizes as a consequence of the edge length constraints).

The defect in the solutions currently provided by the bounded GA is that instead of the generally disc-like expected shape of the convex hull (from squashing the generally spheroid/ovoid shape of the unbounded solutions), an almost cubic convex hull is being formed. This cubic hull is not nearly optimal, as even from a cursory inspection of visualizations of the unbounded and bounded solutions, it can expand a lot more along X and Y.

I surmise this is due to the GA seeing the input as a vector of doubles, and being unable to distinguish between X, Y, and Z coordinates, having no clue of their spatial significance. In creating successive generations, it apparently mixes and matches genes without discrimination. As a result of which, the bounded values from every third gene seem to spread everywhere else.

Answer 1

[Matt J]

This cubic hull is not nearly optimal, as even from a cursory inspection of visualizations of the unbounded and bounded solutions, it can expand a lot more along X and Y.

Well, a test must be made to see whether the fault for that is GA, or whether it is some problem in the constraints or fitness function that you have provided.

It sounds like you could manually alter the solution given by unbounded GA to obtain a solution with a better fitness value. I.e., you could expand it further along X,Y by hand. I would do so, then re-run bounded GA with that hand-crafted solution in the initial population. Then, you can see if GA improves upon that solution or if instead it progresses back toward the solution you don't like. If the latter, you know there is something wrong in the fitness/constraint functions provided by you.

[slaiyer]

penalty = volume(i) / NUM;

for j = 1 : NUM
  if N2(j,3) < bounds(1) || N2(j,3) > bounds(2)
    volume(i) = volume(i) - penalty;
  end
end

I'm currently having to resort to setting a primitive penalty (see above) inside the iteration code for vertices breaching bounds. The solutions are now optimal in all aspects, proving the validity of the code itself (which has already been done for superset intervals in the unbounded version).

However, this approach is hardly "organic", considering that the LB and UB parameters exist exclusively for this purpose. I also tried setting (2 * NUM) linear constraints targeted at the NUM number of Z coordinates, leading to EXACTLY the same behaviour as setting LB and UB. (I even tried setting both the linear constraints and LB/UB together, with all possible valid permutations in the options structure for the sake of completeness.)

Hand-crafting solutions is not an option as this is only part of the big picture that will be completely automated. Although I cannot "accept" your answer as the button so offers, I would like to thank you for your interest in my issue and your efforts towards helping me resolve it.

[Matt J]

Hand-crafting solutions is not an option as this is only part of the big picture that will be completely automated.

Well, it was an option as a test, which is what I proposed it as. Interesting, though, that the penalty approach resolves everything.

[slaiyer]

It doesn't resolve everything, as vertices may still breach bounds, although it's a rare occurrence in the optimized spatial arrangement, and by a minor amount at its worst. This is clearly because the penalty approach is only a discouragement and not a hard limit in any way.

[Matt J]

Well now, you can take that solution and put it in the initial population of bounded GA, like I was suggesting before, and see if it solves the problem rigorously (or at least within TolCon).

Remember, even though it's a "global optimization" algorithm, it can still probably benefit from your help in finding the right region to search.