added regular grid optimization component

docs/source/examples.rst

@@ -25,4 +25,5 @@ You can also `launch them via Binder <https://mybinder.org/v2/git/https%3A%2F%2F
     examples_nblinks/bathymetry_nb
     examples_nblinks/exclusion_zones_nb
     examples_nblinks/mongodb_nb
+    examples_nblinks/regular_grid_optimization_nb
 

docs/source/examples_nblinks/regular_grid_optimization_nb.nblink

+{
+    "path": "../../notebooks/regular_grid_optimization.ipynb"
+}
\ No newline at end of file

examples/docs/example_14_regular_grid_optimization.py

+from py_wake.examples.data.hornsrev1 import V80, Hornsrev1Site
+from py_wake import BastankhahGaussian
+from py_wake.utils.gradients import autograd
+from topfarm import TopFarmProblem
+from topfarm.plotting import XYPlotComp
+from topfarm.easy_drivers import EasyScipyOptimizeDriver
+from topfarm.cost_models.cost_model_wrappers import CostModelComponent
+from topfarm.constraint_components.boundary import XYBoundaryConstraint
+from topfarm.constraint_components.spacing import SpacingConstraint
+
+site = Hornsrev1Site()
+wt = V80()
+D = wt.diameter()
+windFarmModel = BastankhahGaussian(site, wt)
+n_wt = 16
+ec = 0.1
+boundary = [[-800,-50], [1200, -50], [1200,2300], [-800, 2300]]
+
+def aep_fun(x, y):
+    aep = windFarmModel(x, y).aep().sum()
+    return aep
+
+daep = windFarmModel.aep_gradients(gradient_method=autograd, wrt_arg=['x', 'y'])
+aep_comp = CostModelComponent(input_keys=['x', 'y'],
+                              n_wt=n_wt,
+                              cost_function=aep_fun,
+                              cost_gradient_function = daep,
+                              output_keys= ("aep", 0),
+                              output_unit="GWh",
+                              maximize=True,
+                              objective=True)
+
+problem = TopFarmProblem(design_vars={'sx': (3*D, 2*D, 15*D), 
+                                      'sy': (4*D, 2*D, 15*D), 
+                                       'rotation': (50, 0, 90)
+                                      },
+                         constraints=[XYBoundaryConstraint(boundary),
+                                      SpacingConstraint(4*D)],
+                        n_wt = n_wt,
+                        cost_comp=aep_comp,
+                        driver=EasyScipyOptimizeDriver(optimizer='SLSQP', maxiter=200),
+                        plot_comp=XYPlotComp(),
+                        expected_cost=ec,
+                        additional_input={'stagger': 1*D}
+                        )
+
+cost, state, recorder = problem.optimize()
+
+
+

topfarm/_topfarm.py

@@ -33,7 +33,7 @@ from topfarm.recorders import NestedTopFarmListRecorder,\
 from topfarm.mongo_recorder import MongoRecorder
 from topfarm.plotting import NoPlot
 from topfarm.easy_drivers import EasyScipyOptimizeDriver, EasySimpleGADriver, EasyDriverBase
-from topfarm.utils import smart_start
+from topfarm.utils import smart_start, RegularGridComponent
 from topfarm.constraint_components.spacing import SpacingComp
 from topfarm.constraint_components.boundary import BoundaryBaseComp
 from topfarm.constraint_components.penalty_component import PenaltyComponent, PostPenaltyComponent
@@ -81,7 +81,7 @@ class TopFarmProblem(Problem):
                  constraints=[], plot_comp=NoPlot(), record_id=None,
                  expected_cost=1, ext_vars={}, post_constraints=[], approx_totals=False,
                  recorder=None, additional_recorders=None,
-                 n_wt=0):
+                 n_wt=0, additional_input={}):
         """Initialize TopFarmProblem
 
         Parameters
@@ -210,6 +210,8 @@ class TopFarmProblem(Problem):
             self.indeps.add_output(k, v)
         self.ext_vars = ext_vars
 
+        if 'sx' in design_vars:
+            self.model.add_subsystem('regular_grid_comp', RegularGridComponent(design_vars, n_wt, **additional_input), promotes=['*'])
         if cost_comp and isinstance(cost_comp, PostConstraint):
             post_constraints = post_constraints + [cost_comp]
 

topfarm/utils.py

@@ -5,6 +5,8 @@ import multiprocessing
 import threading
 import time
 from tqdm import tqdm
+from openmdao.core.explicitcomponent import ExplicitComponent
+import topfarm
 
 
 def smart_start(XX, YY, ZZ, N_WT, min_space, radius=None, random_pct=0, plot=False, seed=None):
@@ -145,8 +147,164 @@ def smooth_max_gradient(X, alpha, axis=0):
         (1 + alpha * (X - np.expand_dims(smooth_max(X, alpha, axis=axis), axis)))
 
 
+def regular_generic_layout(n_wt, sx, sy, stagger, rotation, ratio=1.0):
+    '''
+    Parameters
+    ----------
+    n_wt : int
+        number of wind turbines
+    sx : float
+        spacing (in turbine diameters or meters) between turbines in x direction
+    sy : float
+        spacing (in turbine diameters or meters) between turbines in y direction
+    stagger : float
+        stagger (in turbine diameters or meters) distance every other turbine column
+    rotation : float
+        rotational angle of the grid in degrees
+    ratio : float
+        ratio between number of columns and number of rows (1.0)
+
+    Returns
+    -------
+    xy : array
+        2D array of x- and y-coordinates (in turbine diameters or meters)
+
+     '''
+
+    if ratio > 1:
+        n_row = int(np.round((n_wt * ratio) ** 0.5))
+        n_col = int(n_wt / n_row)
+    else:
+        n_col = int(np.round((n_wt * ratio) ** 0.5))
+        n_row = int(n_wt / n_col)
+    rest = n_wt - n_col * n_row
+    theta = np.radians(float(rotation))
+    R = np.array([[np.cos(theta), -np.sin(theta)],
+                  [np.sin(theta), np.cos(theta)]])
+    x_grid = np.linspace(0, n_col * float(sx), n_col, endpoint=False)
+    y_grid = np.linspace(0, n_row * float(sy), n_row, endpoint=False)
+    xm, ym = np.meshgrid(x_grid, y_grid)
+    ym[:, 1::2] += stagger
+    x, y = xm.ravel(), ym.ravel()
+    x = np.hstack((x, x_grid[:rest]))
+    y = np.hstack((y, (y[-n_col:-(n_col - rest)] + float(sy))))
+    return np.matmul(R, np.array([x, y]))
+
+
+def regular_generic_layout_gradients(n_wt, sx, sy, stagger, rotation, ratio=1.0):
+    '''
+    Parameters
+    ----------
+    n_wt : int
+        number of wind turbines
+    sx : float
+        spacing (in turbine diameters or meters) between turbines in x direction
+    sy : float
+        spacing (in turbine diameters or meters) between turbines in y direction
+    stagger : float
+        stagger (in turbine diameters or meters) distance every other turbine column
+    rotation : float
+        rotational angle of the grid in degrees
+    ratio : float
+        ratio between number of columns and number of rows (1.0)
+
+    Returns
+    -------
+    dx_dsx, dy_dsx, dx_dsy, dy_dsy, dx_dr, dy_dr : tuple
+        tuple of gradients of x and y with respect to x-spacing, y-spacing and grid rotation
+
+     '''
+
+    if ratio > 1:
+        n_row = int(np.round((n_wt * ratio) ** 0.5))
+        n_col = int(n_wt / n_row)
+    else:
+        n_col = int(np.round((n_wt * ratio) ** 0.5))
+        n_row = int(n_wt / n_col)
+    rest = n_wt - n_col * n_row
+    theta = np.radians(float(rotation))
+    R = np.array([[np.cos(theta), -np.sin(theta)],
+                  [np.sin(theta), np.cos(theta)]])
+
+    x_grid = np.linspace(0, n_col, n_col, endpoint=False)
+    y_grid = np.zeros(n_row)
+    xm, ym = np.meshgrid(x_grid, y_grid)
+    x, y = xm.ravel(), ym.ravel()
+    x = np.hstack((x, x_grid[:rest]))
+    y = np.hstack((y, (y[-n_col:-(n_col - rest)])))
+    dxy_dsx = np.matmul(R, np.array([x, y]))
+    dx_dsx = dxy_dsx[0, :]
+    dy_dsx = dxy_dsx[1, :]
+
+    x_grid = np.zeros(n_col)
+    y_grid = np.linspace(0, n_row, n_row, endpoint=False)
+    xm, ym = np.meshgrid(x_grid, y_grid)
+    x, y = xm.ravel(), ym.ravel()
+    x = np.hstack((x, x_grid[:rest]))
+    y = np.hstack((y, (y[-n_col:-(n_col - rest)] + 1)))
+    dxy_dsy = np.matmul(R, np.array([x, y]))
+    dx_dsy = dxy_dsy[0, :]
+    dy_dsy = dxy_dsy[1, :]
+
+    x_grid = np.linspace(0, n_col * float(sx), n_col, endpoint=False)
+    y_grid = np.linspace(0, n_row * float(sy), n_row, endpoint=False)
+    xm, ym = np.meshgrid(x_grid, y_grid)
+    ym[:, 1::2] += stagger
+    x, y = xm.ravel(), ym.ravel()
+    x = np.hstack((x, x_grid[:rest]))
+    y = np.hstack((y, (y[-n_col:-(n_col - rest)] + float(sy))))
+    dRdr = np.array([[-np.sin(theta), -np.cos(theta)],
+                     [np.cos(theta), -np.sin(theta)]])
+    dx_dr, dy_dr = np.matmul(dRdr, np.array([x, y])) * np.pi / 180
+
+    return dx_dsx, dy_dsx, dx_dsy, dy_dsy, dx_dr, dy_dr
+
+
+class RegularGridComponent(ExplicitComponent):
+    def __init__(self, design_vars, n_wt, **additional_input):
+        super().__init__()
+        self.design_vars = design_vars
+        self.n_wt = n_wt
+        default_vals = {'n_wt': n_wt,
+                        'sx': None,
+                        'sy': None,
+                        'stagger': 0,
+                        'rotation': 0,
+                        'ratio': 1.0}
+        additional_vals = {k: v for k, v in additional_input.items() if k in default_vals}
+        self.default_dict = {**default_vals, **additional_vals}
+
+    def setup(self):
+        for k, v in self.design_vars.items():
+            if k in self.default_dict:
+                self.add_input(k, v[0], units=v[-1])
+                if k in ['sx', 'sy', 'rotation']:
+                    self.declare_partials(topfarm.x_key, k, method='exact')
+                    self.declare_partials(topfarm.y_key, k, method='exact')
+                else:
+                    self.declare_partials(topfarm.x_key, k, method='fd')
+                    self.declare_partials(topfarm.y_key, k, method='fd')
+        self.add_output(topfarm.x_key, val=np.zeros(self.n_wt))
+        self.add_output(topfarm.y_key, val=np.zeros(self.n_wt))
+
+    def compute(self, inputs, outputs):
+        x, y = regular_generic_layout(**{**self.default_dict, **dict(inputs)})
+        outputs[topfarm.x_key] = x
+        outputs[topfarm.y_key] = y
+
+    def compute_partials(self, inputs, J):
+        dx_dsx, dy_dsx, dx_dsy, dy_dsy, dx_dr, dy_dr = regular_generic_layout_gradients(**{**self.default_dict, **dict(inputs)})
+        J[topfarm.x_key, 'sx'] = dx_dsx
+        J[topfarm.x_key, 'sy'] = dx_dsy
+        J[topfarm.y_key, 'sx'] = dy_dsx
+        J[topfarm.y_key, 'sy'] = dy_dsy
+        J[topfarm.x_key, 'rotation'] = dx_dr
+        J[topfarm.y_key, 'rotation'] = dy_dr
+
+
 def main():
     if __name__ == '__main__':
+        plt.close('all')
         N_WT = 30
         min_space = 2.1
 
@@ -156,6 +314,7 @@ def main():
         val = np.sin(XX) + np.sin(YY)
         min_space = 2.1
         xs, ys = smart_start(XX, YY, val, N_WT, min_space)
+        plt.figure()
         c = plt.contourf(XX, YY, val, 100)
         plt.colorbar(c)
         for i in range(N_WT):
@@ -165,5 +324,42 @@ def main():
         plt.axis('equal')
         plt.show()
 
+        plt.figure()
+        xy = regular_generic_layout(n_wt=100, sx=5, sy=4, stagger=2, rotation=35, ratio=1.25)
+
+        dsx = 0.1
+        dxyx = regular_generic_layout(n_wt=100, sx=5 + dsx, sy=4, stagger=2, rotation=35, ratio=1.25)
+        dxy_dsx = (dxyx - xy) / dsx
+        dx_dsx, dy_dsx = dxy_dsx
+
+        dsy = 0.1
+        dxyy = regular_generic_layout(n_wt=100, sx=5, sy=4 + dsy, stagger=2, rotation=35, ratio=1.25)
+        dxy_dsy = (dxyy - xy) / dsy
+        dx_dsy, dy_dsy = dxy_dsy
+
+        dR = 0.1
+        dxyR = regular_generic_layout(n_wt=100, sx=5, sy=4, stagger=2, rotation=35 + dR, ratio=1.25)
+        dxy_dR = (dxyR - xy) / dR
+        dx_dR, dy_dR = dxy_dR
+
+        x, y = xy
+        plt.plot(x, y, '.')
+        plt.axis('equal')
+
+        dx_dsxg, dy_dsxg, dx_dsyg, dy_dsyg, dx_dRg, dy_dRg = regular_generic_layout_gradients(n_wt=100, sx=5, sy=4, stagger=2, rotation=35, ratio=1.25)
+        plt.figure()
+        plt.plot(dx_dsx.ravel(), dx_dsxg.ravel(), '.')
+        plt.figure()
+        plt.plot(dy_dsx.ravel(), dy_dsxg.ravel(), '.')
+        plt.figure()
+        plt.plot(dx_dsy.ravel(), dx_dsyg.ravel(), '.')
+        plt.figure()
+        plt.plot(dy_dsy.ravel(), dy_dsyg.ravel(), '.')
+
+        plt.figure()
+        plt.plot(dx_dR.ravel(), dx_dRg.ravel(), '.')
+        plt.figure()
+        plt.plot(dy_dR.ravel(), dy_dRg.ravel(), '.')
+
 
 main()