Merge branch '68-constraints_as_penalty' into 'master'

Resolve "Implement boundary and spacing as penalty"

Closes #68

See merge request !62

topfarm/_topfarm.py

@@ -186,10 +186,17 @@ class TurbineXYZOptimizationProblem(TopFarmProblem):
         self.min_spacing = min_spacing
 
         spacing_comp = SpacingComp(n_wt, min_spacing)
-        spacing_comp.setup_as_constraints(self)
-
         self.boundary_comp = boundary_comp
-        boundary_comp.setup_as_constraints(self)
+
+        if self.driver.supports['inequality_constraints']:
+            spacing_comp.setup_as_constraints(self)
+            boundary_comp.setup_as_constraints(self)
+            mode = 'fwd'
+        else:
+            spacing_comp.setup_as_penalty(self)
+            boundary_comp.setup_as_penalty(self)
+            mode = 'rev'
+        
 
         do = self.driver.options
         if len(boundary_comp.xy_boundary) > 0:
@@ -231,7 +238,8 @@ class TurbineXYZOptimizationProblem(TopFarmProblem):
             else:
                 plot_comp.n_vertices = 0
 
-        self.setup(check=True, mode='fwd')
+        self.plot_comp = plot_comp
+        self.setup(check=True, mode=mode)
 
     @property
     def turbine_positions(self):

topfarm/constraint_components/boundary_component.py

@@ -58,6 +58,48 @@ class BoundaryBaseComp(ExplicitComponent):
         if len(self.z_boundary):
             problem.model.add_constraint('turbineZ', lower=self.z_boundary[:, 0], upper=self.z_boundary[:, 1])
 
+    def setup_as_penalty(self, problem, penalty=1e20):
+        if len(self.xy_boundary) == 0 and len(self.z_boundary) == 0:
+            return  # no boundary or hub-height constraints
+
+        if len(self.xy_boundary) > 0:
+            subsystem_order = [ss.name for ss in problem.model._static_subsystems_allprocs]
+            problem.model.add_subsystem('xy_bound_comp', self, promotes=['*'])
+            subsystem_order.insert(subsystem_order.index('cost_comp'), 'xy_bound_comp')
+            problem.model.set_order(subsystem_order)
+
+            def xy_boundary_satisfied(inputs):
+                return not np.any(inputs['boundaryDistances'] < self.zeros)
+        else:
+            def xy_boundary_satisfied(inputs):
+                return True
+
+        if len(self.z_boundary):
+            def z_boundary_satisfied(inputs):
+                return not (np.any(inputs['turbineZ'] < self.z_boundary[:, 0]) or np.any(inputs['turbineZ'] > self.z_boundary[:, 1]))
+        else:
+            def z_boundary_satisfied(inputs):
+                return True
+
+        self._cost_comp = problem.cost_comp
+        self._org_setup = self._cost_comp.setup
+        self._org_compute = self._cost_comp.compute
+
+        def new_setup():
+            self._org_setup()
+            if len(self.xy_boundary) > 0:
+                self._cost_comp.add_input('boundaryDistances', val=self.zeros)
+
+        self._cost_comp.setup = new_setup
+
+        def new_compute(inputs, outputs):
+            if xy_boundary_satisfied(inputs) and z_boundary_satisfied(inputs):
+                self._org_compute(inputs, outputs)
+            else:
+                outputs['cost'] = penalty
+        self._cost_comp.compute = new_compute
+        problem._mode = 'rev'
+
 
 class ConvexBoundaryComp(BoundaryBaseComp):
     def __init__(self, n_wt, xy_boundary=None, z_boundary=None, xy_boundary_type='convex_hull'):
@@ -66,6 +108,9 @@ class ConvexBoundaryComp(BoundaryBaseComp):
             self.boundary_type = xy_boundary_type
             self.calculate_boundary_and_normals()
             self.calculate_gradients()
+            self.zeros = np.zeros([self.n_wt, self.nVertices])
+        else:
+            self.zeros = np.zeros([self.n_wt, 0])
 
     def calculate_boundary_and_normals(self):
         if self.boundary_type == 'convex_hull':
@@ -179,7 +224,7 @@ class ConvexBoundaryComp(BoundaryBaseComp):
 
         # Explicitly size output array
         # (vector with positive elements if turbines outside of hull)
-        self.add_output('boundaryDistances', np.zeros([self.n_wt, self.nVertices]),
+        self.add_output('boundaryDistances', self.zeros,
                         desc="signed perpendicular distances from each turbine to each face CCW; + is inside")
 
         self.declare_partials('boundaryDistances', ['turbineX', 'turbineY'])
@@ -227,6 +272,7 @@ class PolygonBoundaryComp(BoundaryBaseComp):
         BoundaryBaseComp.__init__(self, n_wt, xy_boundary=xy_boundary, z_boundary=z_boundary, **kwargs)
 
         self.nTurbines = n_wt
+        self.zeros = np.zeros(self.nTurbines)
         vertices = self.xy_boundary
         self.nVertices = vertices.shape[0]
 
@@ -271,7 +317,7 @@ class PolygonBoundaryComp(BoundaryBaseComp):
 
         # Explicitly size output array
         # (vector with positive elements if turbines outside of hull)
-        self.add_output('boundaryDistances', np.zeros(self.nTurbines),
+        self.add_output('boundaryDistances', self.zeros,
                         desc="signed perpendicular distances from each turbine to each face CCW; + is inside")
 
         self.declare_partials('boundaryDistances', ['turbineX', 'turbineY'])

topfarm/constraint_components/spacing_component.py

@@ -19,6 +19,32 @@ class SpacingComp(ExplicitComponent):
             zero = np.zeros(int(((self.n_wt - 1.) * self.n_wt / 2.)))
             problem.model.add_constraint('wtSeparationSquared', lower=zero + (self.min_spacing)**2)
 
+    def setup_as_penalty(self, problem, penalty=1e20):
+        if self.min_spacing is not None:
+            subsystem_order = [ss.name for ss in problem.model._static_subsystems_allprocs]
+            problem.model.add_subsystem('spacing_comp', self, promotes=['*'])
+            subsystem_order.insert(subsystem_order.index('cost_comp'), 'spacing_comp')
+            problem.model.set_order(subsystem_order)
+            
+            zero = np.zeros(int(((self.n_wt - 1.) * self.n_wt / 2.)))
+
+            self._cost_comp = problem.cost_comp
+            self._org_setup = self._cost_comp.setup
+            self._org_compute = self._cost_comp.compute
+
+            def new_setup():
+                self._org_setup()
+                self._cost_comp.add_input('wtSeparationSquared', val=zero)
+
+            self._cost_comp.setup = new_setup
+
+            def new_compute(inputs, outputs):
+                if np.any(inputs['wtSeparationSquared'] < zero + self.min_spacing**2):
+                    outputs['cost'] = penalty
+                else:
+                    self._org_compute(inputs, outputs)
+            self._cost_comp.compute = new_compute
+
     def setup(self):
 
         # set finite difference options (fd used for testing only)

Colonel @ 315a896d

-Subproject commit db2087ac3f1075028ffd5c0b5de5732b6664086a
+Subproject commit 315a896dfe20475ca1524523a021e7dedb3e64cd

topfarm/tests/test_constraint/test_BoundaryBaseComp.py

@@ -6,62 +6,96 @@ from topfarm import TopFarm
 import pytest
 from topfarm._topfarm import TurbineXYZOptimizationProblem
 from topfarm.constraint_components.boundary_component import BoundaryComp
+from topfarm.plotting import NoPlot
+from openmdao.drivers.genetic_algorithm_driver import SimpleGADriver
+import warnings
 
+initial = np.array([[0, 0, 0], [6, 0, 0], [6, -10, 0]])  # initial turbine layouts
+optimal = np.array([[2.5, -3], [6, -7], [4.5, -3]])  # desired turbine layouts
+boundary = np.array([(0, 0), (6, 0), (6, -10), (0, -10)])  # turbine boundaries
+desired = np.array([[3, -3, 1], [7, -7, 2], [4, -3, 3]])  # desired turbine layouts
 
-def testSquare():
-    optimal = [(0, 0)]
-    boundary = [(0, 0), (1, 3)]
-    tf = TopFarm(optimal, DummyCost(optimal), 2, boundary=boundary, boundary_type='square', record_id=None)
-    np.testing.assert_array_equal(tf.xy_boundary, [[-1, 0],
-                                                   [2, 0],
-                                                   [2, 3],
-                                                   [-1, 3],
-                                                   [-1, 0]])
-
-
-def testRectangle():
-    optimal = [(0, 0)]
-    boundary = [(0, 0), (1, 3)]
-    tf = TopFarm(optimal, DummyCost(optimal), 2, boundary=boundary, boundary_type='rectangle', record_id=None)
-    np.testing.assert_array_equal(tf.xy_boundary, [[0, 0],
-                                                   [1, 0],
-                                                   [1, 3],
-                                                   [0, 3],
-                                                   [0, 0]])
-
-
-def testConvexHull():
-    optimal = [(0, 0)]
-    boundary = [(0, 0), (1, 1), (2, 0), (2, 2), (0, 2)]
-    tf = TopFarm(optimal, DummyCost(optimal), 2, boundary=boundary, boundary_type='convex_hull', record_id=None)
-    np.testing.assert_array_equal(tf.xy_boundary, [[0, 0],
-                                                   [2, 0],
-                                                   [2, 2],
-                                                   [0, 2],
-                                                   [0, 0]])
-
-
-def testNotImplemented():
-    optimal = [(0, 0)]
-    boundary = [(0, 0), (1, 1), (2, 0), (2, 2), (0, 2)]
-    with pytest.raises(NotImplementedError, match="Boundary type 'Something' is not implemented"):
-        TopFarm(optimal, DummyCost(optimal), 2, boundary=boundary, boundary_type='Something', record_id=None)
-
-
-def test_z_boundary():
-    optimal = [(0, 0)]
-    tf = TurbineXYZOptimizationProblem(DummyCost(optimal), optimal, BoundaryComp(2, None, [70, 90]))
-    np.testing.assert_array_equal(tf.z_boundary, [[70, 90],
-                                                  [70, 90]])
-
-
-def test_xyz_boundary():
-    optimal = [(0, 0)]
-    boundary = [(0, 0), (1, 3)]
-    tf = TurbineXYZOptimizationProblem(DummyCost(optimal), optimal, BoundaryComp(1, boundary, [70, 90], xy_boundary_type='rectangle'))
-    np.testing.assert_array_equal(tf.xy_boundary, [[0, 0],
-                                                   [1, 0],
-                                                   [1, 3],
-                                                   [0, 3],
-                                                   [0, 0]])
-    np.testing.assert_array_equal(tf.z_boundary, [[70, 90]])
+
+def test_setup_as_constraint_xy():
+    from topfarm.cost_models.dummy import DummyCostPlotComp
+
+    #plot_comp = DummyCostPlotComp(desired)
+    plot_comp = NoPlot()
+
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired[:, :2], ['turbineX', 'turbineY']), initial[:, :2],
+                                       boundary_comp=BoundaryComp(len(initial), boundary, None),
+                                       plot_comp=plot_comp)
+
+    tf.optimize()
+    tb_pos = tf.turbine_positions[:, :2]
+    tf.plot_comp.show()
+    tol = 1e-4
+    assert tb_pos[1][0] < 6 + tol  # check within border
+
+
+def test_setup_as_constraint_z():
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired[:, 2:], ['turbineZ']), initial,
+                                       boundary_comp=BoundaryComp(len(initial), None, [0, 2]))
+    tf.optimize()
+    assert np.all(tf.turbine_positions[:, 2]) <= 2  # check within height limit
+
+
+def test_setup_as_constraint_xyz():
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired, ['turbineX', 'turbineY', 'turbineZ']), initial,
+                                       boundary_comp=BoundaryComp(len(initial), boundary, [0, 2]))
+    tf.optimize()
+    tb_pos = tf.turbine_positions
+    tol = 1e-4
+    assert tb_pos[1][0] < 6 + tol  # check within border
+    assert np.all(tf.turbine_positions[:, 2]) <= 2  # check within height limit
+
+
+def test_setup_as_penalty_xy():
+    driver = SimpleGADriver()
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired, ['turbineX', 'turbineY']), initial,
+                                       boundary_comp=BoundaryComp(len(initial), boundary, None),
+                                       driver=driver)
+
+    # check normal result if spacing constraint is satisfied
+    assert tf.evaluate()[0] == 121
+    # check penalized result if spacing constraint is not satisfied
+    assert tf.evaluate({'turbineX': [2.5, 7, 4.5], 'turbineY': [-3., -7., -3.], 'turbineZ': [0., 0., 0.]})[0] == 1e20
+
+
+def test_setup_as_penalty_z():
+    driver = SimpleGADriver()
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired[:,2:], ['turbineZ']), initial,
+                                       boundary_comp=BoundaryComp(3, None, [0, 2]),
+                                       driver=driver)
+
+    # check normal result if spacing constraint is satisfied
+    assert tf.evaluate()[0] == 14
+    # check penalized result if spacing constraint is not satisfied
+    assert tf.evaluate({'turbineX': [2.5, 7, 4.5], 'turbineY': [-3., -7., -3.], 'turbineZ': [0., 0., 3.]})[0] == 1e20
+    assert tf.evaluate({'turbineX': [2.5, 7, 4.5], 'turbineY': [-3., -7., -3.], 'turbineZ': [-1., 0., 0.]})[0] == 1e20
+
+
+def test_setup_as_penalty_xyz():
+    driver = SimpleGADriver()
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired, ['turbineX', 'turbineY', 'turbineZ']), initial,
+                                       boundary_comp=BoundaryComp(len(initial), boundary, [0, 2]),
+                                       driver=driver)
+
+    # check normal result if spacing constraint is satisfied
+    assert tf.evaluate()[0] == 135
+    # check penalized result if spacing constraint is not satisfied
+    assert tf.evaluate({'turbineX': [2.5, 7, 4.5], 'turbineY': [-3., -7., -3.], 'turbineZ': [0., 0., 0.]})[0] == 1e20  # outside border
+    assert tf.evaluate({'turbineX': [2.5, 7, 4.5], 'turbineY': [-3., -7., -3.], 'turbineZ': [0., 0., 3.]})[0] == 1e20  # above limit
+    assert tf.evaluate({'turbineX': [2.5, 7, 4.5], 'turbineY': [-3., -7., -3.], 'turbineZ': [-1., 0., 0.]})[0] == 1e20  # below limit
+
+def test_setup_as_penalty_none():
+    driver = SimpleGADriver()
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired, ['turbineX', 'turbineY']), initial,
+                                       boundary_comp=BoundaryComp(len(initial), None, None),
+                                       driver=driver)
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        # check normal result if spacing constraint is satisfied
+        assert tf.evaluate()[0] == 121
+        # check penalized result if spacing constraint is not satisfied
+        assert tf.evaluate({'turbineX': [2.5, 7, 4.5], 'turbineY': [-3., -7., -3.], 'turbineZ': [0., 0., 0.]})[0] == .5

topfarm/tests/test_constraint/test_spacingComp.py

+from topfarm import TopFarm
+import numpy as np
+import pytest
+from topfarm.cost_models.dummy import DummyCost
+from topfarm.plotting import NoPlot
+from topfarm.easy_drivers import EasyScipyOptimizeDriver, EasyPyOptSparseIPOPT
+from topfarm.tests import npt, uta
+from topfarm._topfarm import TurbineXYZOptimizationProblem
+from topfarm.constraint_components.boundary_component import BoundaryComp
+from openmdao.drivers.genetic_algorithm_driver import SimpleGADriver
+import warnings
+
+
+initial = np.array([[6, 0], [6, -8], [1, 1]])  # initial turbine layouts
+optimal = np.array([[2.5, -3], [6, -7], [4.5, -3]])  # desired turbine layouts
+boundary = np.array([(0, 0), (6, 0), (6, -10), (0, -10)])  # turbine boundaries
+desired = np.array([[3, -3], [7, -7], [4, -3]])  # desired turbine layouts
+
+
+def test_spacing():
+    from topfarm.cost_models.dummy import DummyCostPlotComp
+
+    #plot_comp = DummyCostPlotComp(desired)
+    plot_comp = NoPlot()
+
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired, ['turbineX', 'turbineY']), initial,
+                                       boundary_comp=BoundaryComp(len(initial), boundary, None),
+                                       min_spacing=2, plot_comp=plot_comp)
+
+    tf.evaluate()
+    tf.optimize()
+    tb_pos = tf.turbine_positions[:, :2]
+    tf.plot_comp.show()
+    tol = 1e-4
+    assert sum((tb_pos[2] - tb_pos[0])**2) > 2**2 - tol  # check min spacing
+
+
+def test_spacing_as_penalty():
+
+    driver = SimpleGADriver()
+
+    tf = TurbineXYZOptimizationProblem(DummyCost(desired, ['turbineX', 'turbineY']), initial,
+                                       boundary_comp=BoundaryComp(len(initial), None, None),
+                                       min_spacing=2, driver=driver)
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        # check normal result if spacing constraint is satisfied
+        assert tf.evaluate()[0] == 45
+        # check penalized result if spacing constraint is not satisfied
+        assert tf.evaluate({'turbineX': [3, 7, 4.], 'turbineY': [-3., -7., -3.], 'turbineZ': [0., 0., 0.]})[0] == 1e20

topfarm/tests/test_topfarm_problems/test_InitialXYZOptimizationProblem.py

@@ -13,14 +13,14 @@ from topfarm.constraint_components.boundary_component import BoundaryComp
 @pytest.fixture
 def get_tf():
     def get_InitialXYZOptimizationProblem(driver,
-                                          cost_comp=DummyCost([(1, 0, 4),
-                                                               (0, 1, 3)]),
                                           min_spacing=None,
                                           turbineXYZ=[[0, 0, 0],
                                                       [2, 2, 2]],
                                           xy_boundary=[(10, 6), (11, 8)],
                                           xy_boundary_type='rectangle',
                                           z_boundary=[3, 4]):
+        cost_comp = DummyCost([(1, 0, 4),
+                               (0, 1, 3)])
         return InitialXYZOptimizationProblem(
             cost_comp, turbineXYZ,
             BoundaryComp(len(turbineXYZ), xy_boundary, z_boundary, xy_boundary_type),

topfarm/tests/test_topfarm_problems/test_TurbineTypeOptimizationProblem.py

+from topfarm._topfarm import TurbineTypeOptimizationProblem
+from openmdao.drivers.doe_generators import FullFactorialGenerator
+from openmdao.drivers.doe_driver import DOEDriver
+from topfarm.cost_models.dummy import DummyCost
+from topfarm.tests import npt
+
+optimal = [[1],[0]]
+
+
+def test_turbineType_optimization():
+    cost_comp = DummyCost(
+        optimal_state=optimal,
+        inputs=['turbineType'])
+    tf = TurbineTypeOptimizationProblem(
+        cost_comp=cost_comp,
+        turbineTypes=[0, 0], lower=0, upper=1,
+        driver=DOEDriver(FullFactorialGenerator(2)))
+    cost, state, _ = tf.optimize()
+    assert cost == 0
+    npt.assert_array_equal(state['turbineType'], [1, 0])