first attempt to support polygon

topfarm/constraint_components/boundary_component.py

 import numpy as np
 from openmdao.api import Group, IndepVarComp, ExecComp, ExplicitComponent, Problem
 from scipy.spatial import ConvexHull
+from shapely.geometry.polygon import Polygon, LinearRing
+from shapely.geometry.point import Point
+from shapely.geometry.multipoint import MultiPoint
 
 
 class BoundaryComp(ExplicitComponent):
@@ -14,7 +17,6 @@ class BoundaryComp(ExplicitComponent):
         self.nVertices = self.vertices.shape[0]
 
     def calculate_boundary_and_normals(self, vertices, boundary_type):
-
         if boundary_type == 'convex_hull':
             # find the points that actually comprise a convex hull
             hull = ConvexHull(list(vertices))
@@ -147,3 +149,56 @@ class BoundaryComp(ExplicitComponent):
         # return Jacobian dict
         partials['boundaryDistances', 'turbineX'] = dfaceDistance_dx
         partials['boundaryDistances', 'turbineY'] = dfaceDistance_dy
+
+
+
+class PolygonBoundaryComp(BoundaryComp):
+
+    def __init__(self, vertices, nTurbines):
+
+        super(BoundaryComp, self).__init__()
+
+        self.nTurbines = nTurbines
+        self.vertices = np.array(vertices)
+        self.nVertices = self.vertices.shape[0]
+        self.polygon = Polygon(self.vertices)
+        self.linearRing = LinearRing(self.polygon.exterior.coords)
+
+
+    def setup(self):
+
+        # Explicitly size input arrays
+        self.add_input('turbineX', np.zeros(self.nTurbines), units='m',
+                       desc='x coordinates of turbines in global ref. frame')
+        self.add_input('turbineY', np.zeros(self.nTurbines), units='m',
+                       desc='y coordinates of turbines in global ref. frame')
+
+        # Explicitly size output array
+        # (vector with positive elements if turbines outside of hull)
+        self.add_output('boundaryDistances', np.zeros([self.nTurbines]),
+                        desc="signed shortest distance to boundary; + is inside")
+
+        # self.declare_partials('boundaryDistances', ['turbineX', 'turbineY'])
+        self.declare_partials('boundaryDistances', ['turbineX', 'turbineY'], method='fd')
+
+  
+    def compute(self, inputs, outputs):
+
+        turbineX = inputs['turbineX']
+        turbineY = inputs['turbineY']
+
+        distances = []
+        for x, y in zip(turbineX, turbineY):
+            point = Point(x, y)
+            d = self.linearRing.project(point)
+            p = self.linearRing.interpolate(d)
+            closest_point_on_boundary = list(p.coords)[0]
+            distances.append([-1, 1][self.polygon.contains(point)] * point.distance(p))
+        print (np.array([turbineX,turbineY]).T)
+        print (distances)
+        print ()
+        outputs['boundaryDistances'] = distances
+        
+    def compute_partials(self, inputs, partials):
+        return
+  
\ No newline at end of file

topfarm/plotting.py

@@ -18,8 +18,8 @@ def mypause(interval):
                 canvas.draw()
             canvas.start_event_loop(interval)
             return
-        
-        
+
+
 class PlotComp(ExplicitComponent):
     """
     Evaluates the equation f(x,y) = (x-3)^2 + xy + (y+4)^2 - 3.
@@ -64,7 +64,7 @@ class PlotComp(ExplicitComponent):
         y = inputs['turbineY']
         cost = inputs['cost']
         self.history = [(x.copy(), y.copy())] + self.history[:self.memory]
-        
+
         boundary = inputs['boundary']
         self.init_plot(boundary)
         plt.title(cost)
@@ -74,9 +74,9 @@ class PlotComp(ExplicitComponent):
             plt.plot(history_arr[:, 0, i], history_arr[:, 1, i], '.-', color=c, lw=1)
             plt.plot(x_, y_, 'o', color=c, ms=5)
             plt.plot(x_, y_, 'x' + 'k', ms=4)
-        
-        if self.counter ==0:
+
+        if self.counter == 0:
             plt.pause(.01)
         mypause(0.01)
-        
+
         self.counter += 1

topfarm/topfarm.py

@@ -3,7 +3,8 @@ import time
 from openmdao.api import Problem, ScipyOptimizeDriver, IndepVarComp
 
 import numpy as np
-from topfarm.constraint_components.boundary_component import BoundaryComp
+from topfarm.constraint_components.boundary_component import BoundaryComp,\
+    PolygonBoundaryComp
 from topfarm.constraint_components.spacing_component import SpacingComp
 import warnings
 
@@ -20,9 +21,12 @@ class TopFarm(object):
                  driver_options={'optimizer': 'SLSQP'}):
 
         self.initial_positions = turbines = np.array(turbines)
-        
+
         n_wt = turbines.shape[0]
-        self.boundardy_comp = BoundaryComp(boundary, n_wt, boundary_type)
+        if boundary_type == 'polygon':
+            self.boundardy_comp = PolygonBoundaryComp(boundary, n_wt)
+        else:
+            self.boundardy_comp = BoundaryComp(boundary, n_wt, boundary_type)
         self.problem = prob = Problem()
         indeps = prob.model.add_subsystem('indeps', IndepVarComp(), promotes=['*'])
         indeps.add_output('turbineX', turbines[:, 0], units='m')
@@ -45,6 +49,7 @@ class TopFarm(object):
             plot_comp.n_wt = n_wt
             plot_comp.n_vertices = self.boundardy_comp.vertices.shape[0]
             prob.model.add_subsystem('plot_comp', plot_comp, promotes=['*'])
+
         self.plot_comp = plot_comp
         prob.model.add_constraint('wtSeparationSquared', lower=np.zeros(int(((n_wt - 1.) * n_wt / 2.))) + (min_spacing)**2)
         prob.model.add_constraint('boundaryDistances', lower=np.zeros(self.boundardy_comp.nVertices * n_wt))
@@ -71,13 +76,13 @@ class TopFarm(object):
         with warnings.catch_warnings():  # suppress OpenMDAO/SLSQP warnings
             warnings.filterwarnings('ignore', "Inefficient choice of derivative mode.  You chose 'rev' for a problem with")
             self.problem.run_model()
-        print ("Evaluated in\t%.3fs" % (time.time() - t))
+        print("Evaluated in\t%.3fs" % (time.time() - t))
         return self.get_cost(), self.turbine_positions
 
     def optimize(self):
         t = time.time()
         self.problem.run_driver()
-        print ("Optimized in\t%.3fs" % (time.time() - t))
+        print("Optimized in\t%.3fs" % (time.time() - t))
         return np.array([self.problem['turbineX'], self.problem['turbineY']]).T
 
     def get_cost(self):
@@ -86,7 +91,7 @@ class TopFarm(object):
     @property
     def boundary(self):
         b = self.boundardy_comp.vertices
-        return np.r_[b,b[:1]]
+        return np.r_[b, b[:1]]
 
     @property
     def turbine_positions(self):
@@ -96,7 +101,7 @@ class TopFarm(object):
 if __name__ == '__main__':
     from topfarm.cost_models.dummy import DummyCostPlotComp, DummyCost
     from topfarm.plotting import PlotComp
-    
+
     n_wt = 4
     random_offset = 5
     optimal = [(3, -3), (7, -7), (4, -3), (3, -7), (-3, -3), (-7, -7), (-4, -3), (-3, -7)][:n_wt]