diff --git a/topfarm/_topfarm.py b/topfarm/_topfarm.py index 92485e04a31cd9a5803b0e837edcedb82de84ea3..a224149523dd56749d58bd9f669cc0c407d7cf80 100644 --- a/topfarm/_topfarm.py +++ b/topfarm/_topfarm.py @@ -3,6 +3,7 @@ from topfarm.constraint_components.boundary_component import BoundaryComp,\ from topfarm.constraint_components.spacing_component import SpacingComp from topfarm.plotting import PlotComp from topfarm.utils import pos_from_case, latest_id +from topfarm.utils import shuffle_positions as spos import os import time import numpy as np @@ -12,6 +13,7 @@ with warnings.catch_warnings(): from openmdao.api import Problem, ScipyOptimizeDriver, IndepVarComp, \ SqliteRecorder + class TopFarm(object): """Optimize wind farm layout in terms of - Position of turbines @@ -24,6 +26,7 @@ class TopFarm(object): boundary_type='convex_hull', plot_comp=None, driver=ScipyOptimizeDriver(), record=False, case_recorder_dir=os.getcwd(), rerun_case_id=None): + self.min_spacing = min_spacing if rerun_case_id is None: self.initial_positions = turbines = np.array(turbines) elif rerun_case_id is 'latest': @@ -34,6 +37,7 @@ class TopFarm(object): else: self.initial_positions = turbines = pos_from_case(rerun_case_id) n_wt = turbines.shape[0] + self.n_wt = n_wt if boundary_type == 'polygon': self.boundary_comp = PolygonBoundaryComp(boundary, n_wt) else: @@ -44,13 +48,15 @@ class TopFarm(object): min_x, min_y = self.boundary_comp.vertices.min(0) mean_x, mean_y = self.boundary_comp.vertices.mean(0) design_var_kwargs = {} - - if 'optimizer' in driver.options and driver.options['optimizer'] == 'SLSQP': + do = driver.options + if 'optimizer' in do and do['optimizer'] == 'SLSQP': min_x, min_y, mean_x, mean_y = 0, 0, 1, 1 # scaling disturbs SLSQP # Default +/- sys.float_info.max does not work for SLSQP design_var_kwargs = {'lower': np.nan, 'upper': np.nan} - indeps.add_output('turbineX', turbines[:, 0], units='m', ref0=min_x, ref=mean_x) - indeps.add_output('turbineY', turbines[:, 1], units='m', ref0=min_y, ref=mean_y) + indeps.add_output('turbineX', turbines[:, 0], units='m', ref0=min_x, + ref=mean_x) + indeps.add_output('turbineY', turbines[:, 1], units='m', ref0=min_y, + ref=mean_y) indeps.add_output('boundary', self.boundary_comp.vertices, units='m') prob.model.add_subsystem('cost_comp', cost_comp, promotes=['*']) prob.driver = driver @@ -111,7 +117,8 @@ class TopFarm(object): def evaluate_gradients(self): t = time.time() - res = self.problem.compute_totals(['cost'], wrt=['turbineX', 'turbineY'], return_format='dict') + res = self.problem.compute_totals(['cost'], wrt=['turbineX', + 'turbineY'], return_format='dict') print("Gradients evaluated in\t%.3fs" % (time.time() - t)) return res @@ -133,33 +140,52 @@ class TopFarm(object): def turbine_positions(self): return np.array([self.problem['turbineX'], self.problem['turbineY']]).T - - def post_process(self, anim_time=10, verbose=True): - if self.plot_comp.animate: - self.plot_comp.run_animate(anim_time, verbose) + def clean(self): for file in os.listdir(self.plot_comp.temp): if file.startswith('plot_') and file.endswith('.png'): - os.remove(os.path.join(self.plot_comp.temp,file)) + os.remove(os.path.join(self.plot_comp.temp, file)) + + def shuffle_positions(self, shuffle_type='rel', n_iter=1000, + step_size=0.1, pad=1.1, offset=5, plot=False, + verbose=False): + if shuffle_type is not None: + turbines = spos(self.boundary, self.n_wt, self.min_spacing, + self.turbine_positions, shuffle_type, n_iter, + step_size, pad, offset, plot, verbose) + self.problem['turbineX'] = turbines.T[0] + self.problem['turbineY'] = turbines.T[1] + + def animate(self, anim_time=10, verbose=False): + if self.plot_comp.animate: + self.plot_comp.run_animate(anim_time, verbose) + else: + if verbose: + print('Animation requested but was not enabled for this ' + 'optimization. Set plot_comp.animate = True to enable') def try_me(): if __name__ == '__main__': from topfarm.cost_models.dummy import DummyCostPlotComp, DummyCost - n_wt = 4 + n_wt = 20 random_offset = 5 - optimal = [(3, -3), (7, -7), (4, -3), (3, -7), (-3, -3), (-7, -7), (-4, -3), (-3, -7)][:n_wt] + optimal = [(3, -3), (7, -7), (4, -3), (3, -7), (-3, -3), (-7, -7), + (-4, -3), (-3, -7)][:n_wt] rotorDiameter = 1.0 minSpacing = 2.0 - turbines = np.array(optimal) + np.random.randint(-random_offset, random_offset, (n_wt, 2)) plot_comp = DummyCostPlotComp(optimal) - plot_comp.animate = True +# plot_comp.animate = True boundary = [(0, 0), (6, 0), (6, -10), (0, -10)] - tf = TopFarm(turbines, DummyCost(optimal), minSpacing * rotorDiameter, boundary=boundary, plot_comp=plot_comp) + tf = TopFarm(optimal, DummyCost(optimal), minSpacing * rotorDiameter, + boundary=boundary, plot_comp=plot_comp) # tf.check() + tf.shuffle_positions(shuffle_type='abs', offset=random_offset) tf.optimize() - tf.post_process() + tf.animate() + tf.clean() + try_me() diff --git a/topfarm/tests/topfarm/test_utils.py b/topfarm/tests/topfarm/test_utils.py index c4c8400ac2d252d3a3abc78ac9e0dbec3a0c4fff..fb0339a5f1818211c7b5cf38033465dc77f3d44f 100644 --- a/topfarm/tests/topfarm/test_utils.py +++ b/topfarm/tests/topfarm/test_utils.py @@ -5,7 +5,7 @@ from topfarm.cost_models.cost_model_wrappers import CostModelComponent import pytest import os -from topfarm.utils import pos_from_case, latest_id, _random_positions +from topfarm.utils import pos_from_case, latest_id, _shuffle_positions_abs thisdir = os.path.dirname(os.path.abspath(__file__)) turbines = np.array([[ 2.4999377 , -2.99987763], @@ -58,13 +58,13 @@ def testlatest_id(): ref_path = os.path.join(path,'cases_20180621_111710.sql') assert latest_id(path) == ref_path -def test_random_positions(): - turbines2 = _random_positions(x, y, boundary, n_wt, n_iter, step_size, +def test_shuffle_positions_abs(): + turbines2 = _shuffle_positions_abs(x, y, boundary, n_wt, n_iter, step_size, min_space, pad, plot, verbose) np.testing.assert_allclose(turbines2, turbines2_ref) if __name__ == '__main__': # testpos_from_case() # testlatest_id() -# test_random_positions() +# test_shuffle_positions_abs() pass \ No newline at end of file diff --git a/topfarm/utils.py b/topfarm/utils.py index 18894baf5df347bc8b0a6adb317516c95a7513a5..70a8e38b53b164ec8c853553e6ea4d3b20060d03 100644 --- a/topfarm/utils.py +++ b/topfarm/utils.py @@ -2,7 +2,6 @@ from topfarm.constraint_components.boundary_component import \ PolygonBoundaryComp from topfarm.constraint_components.spacing_component import SpacingComp import os -from copy import copy import numpy as np from openmdao.api import CaseReader import matplotlib.pyplot as plt @@ -46,51 +45,76 @@ def latest_id(case_recorder_dir): return latest -def random_positions(boundary, n_wt, n_iter, step_size, min_space, - pad=1.1, plot=False, verbose=True): +def shuffle_positions(boundary, n_wt, min_space, init_pos, shuffle_type='abs', + n_iter=1000, step_size=0.1, pad=1.1, offset=0.5, + plot=True, verbose=True): ''' Input: boundary: list of tuples, e.g.: [(0, 0), (6, 1), (7, -11), (-1, -10)] n_wt: number of wind turbines + min_space: the minimum spacing between turbines + init_pos: inital positions that suffeling is based off of if the + shuffle_type requires it + shuffle_type: + 'abs': absolute random positions that respect the boundary, + and aims to respect the minimum spacing + 'rel': moves each turbine with an offset compared to the + initial positions - not implemented yet. n_iter: number of iterations allowed to try and satisfy the minimum spacing constraint step_size: the multiplier on the spacing gradient that the turbines are moved in each step - min_space: the minimum spacing between turbines pad: the multiplier on the boundary gradient plot: plot the generated random layout verbose: print helpful text to the console Returns an array of xy coordinates of the wind turbines ''' - x, y = map(_random, np.array(boundary).T) - turbines = _random_positions(x, y, boundary, n_wt, n_iter, step_size, - min_space, pad, plot, verbose) + if shuffle_type == 'abs': + def _random(b): + return np.random.rand(n_wt) * (max(b) - min(b)) + min(b) + x, y = map(_random, np.array(boundary).T) + turbines = _shuffle_positions_abs(x, y, boundary, n_wt, n_iter, + step_size, min_space, pad, plot, + verbose) + elif shuffle_type == 'rel': + turbines = _shuffle_postions_rel(init_pos, offset, boundary, n_wt, + plot) + return turbines + + +def _shuffle_postions_rel(init_pos, offset, boundary, n_wt, plot): + turbineX = init_pos[:, 0] + turbineY = init_pos[:, 1] + ba = np.array(boundary).T + turbines = np.array(init_pos) + np.random.rand(n_wt, 2)*2*offset-offset + if plot: + plt.cla() + plt.plot(ba[0], ba[1]) + plt.plot(turbineX, turbineY, '.') + plt.plot(turbines.T[0], turbines.T[1], 'o') return turbines -def _random_positions(x, y, boundary, n_wt, n_iter, step_size, min_space, - pad, plot, verbose): +def _shuffle_positions_abs(turbineX, turbineY, boundary, n_wt, n_iter, + step_size, min_space, pad, plot, verbose): boundary_comp = PolygonBoundaryComp(boundary, n_wt) spacing_comp = SpacingComp(nTurbines=n_wt) - turbineX = copy(x) - turbineY = copy(y) min_space2 = min_space**2 ba = np.array(boundary).T - bx = np.append(ba[0], ba[0][0]) - by = np.append(ba[1], ba[1][0]) if plot: plt.cla() - plt.plot(bx, by) - plt.plot(x, y, '.') + plt.plot(ba[0], ba[1]) + plt.plot(turbineX, turbineY, '.') for j in range(n_iter): dist = spacing_comp._compute(turbineX, turbineY) dx, dy = spacing_comp._compute_partials(turbineX, turbineY) index = np.argmin(dist) - if dist[index] < min_space2: + if dist[index] < min_space2 or j == 0: turbineX += dx[index]*step_size turbineY += dy[index]*step_size - turbineX, turbineY = _contain(n_wt, turbineX, turbineY, - boundary_comp, pad) + turbineX, turbineY = _move_inside_boundary(n_wt, turbineX, + turbineY, boundary_comp, + pad) else: if verbose: print('Obtained required spacing after {} iterations'.format( @@ -107,11 +131,7 @@ def _random_positions(x, y, boundary, n_wt, n_iter, step_size, min_space, return turbines -def _random(b): - return np.random.rand(n_wt) * (max(b) - min(b)) + min(b) - - -def _contain(n_wt, turbineX, turbineY, boundary_comp, pad): +def _move_inside_boundary(n_wt, turbineX, turbineY, boundary_comp, pad): for i in range(0, n_wt): dng = boundary_comp.calc_distance_and_gradients(turbineX, turbineY) dist = dng[0][i] @@ -135,14 +155,11 @@ if __name__ == '__main__': latest_id = latest_id(case_recorder_dir) print(latest_id) - boundary = [(0, 0), (6, 1), (7, -11), (-1, -10)] + boundary = [(0, 0), (6, 1), (7, -11), (-1, -10), (0, 0)] n_wt = 20 - n_iter = 1000 - step_size = 0.1 + init_pos = np.column_stack((np.random.randint(0, 6, (n_wt)), + np.random.randint(-10, 0, (n_wt)))) min_space = 2.1 - pad = 1.01 - plot = True - verbose = True - turbines = random_positions(boundary, n_wt, n_iter, step_size, min_space, - pad, plot, verbose) + turbines = shuffle_positions(boundary, n_wt, min_space, init_pos, + shuffle_type='rel') print(turbines)