diff --git a/topfarm/constraint_components/boundary_component.py b/topfarm/constraint_components/boundary_component.py
index 7548bf427c0e878a7f0e7bbfab33c52a67b5c381..83e51da88f69b46e70e926b451a22e9082e75c65 100644
--- a/topfarm/constraint_components/boundary_component.py
+++ b/topfarm/constraint_components/boundary_component.py
@@ -1,6 +1,9 @@
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
diff --git a/topfarm/plotting.py b/topfarm/plotting.py
index bfee8ae4ad47cf8c5c70d2c70c3020cb7f1cf4ec..5dac3a8246ebdcd74697672f1eb2046e0ba80a47 100644
--- a/topfarm/plotting.py
+++ b/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
diff --git a/topfarm/topfarm.py b/topfarm/topfarm.py
index e019117d044882de85ed33078a5168ba0a1ceac1..668262870dc962eaa319cce5abd958aed8ad37c5 100644
--- a/topfarm/topfarm.py
+++ b/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]