diff --git a/topfarm/constraint_components/spacing_component.py b/topfarm/constraint_components/spacing_component.py index ec97351d85f0706322a4874abbcb631996184596..7d05698157636ec8121f3d903aa72ba924e07159 100644 --- a/topfarm/constraint_components/spacing_component.py +++ b/topfarm/constraint_components/spacing_component.py @@ -33,34 +33,39 @@ class SpacingComp(ExplicitComponent): self.declare_partials('wtSeparationSquared', ['turbineX', 'turbineY']) def compute(self, inputs, outputs): - # print 'in dist const' - turbineX = inputs['turbineX'] turbineY = inputs['turbineY'] + separation_squared = self._compute(turbineX, turbineY) + outputs['wtSeparationSquared'] = separation_squared + + + def _compute(self, turbineX, turbineY): nTurbines = self.nTurbines separation_squared = np.zeros(int((nTurbines - 1) * nTurbines / 2)) - k = 0 for i in range(0, nTurbines): for j in range(i + 1, nTurbines): separation_squared[k] = (turbineX[j] - turbineX[i])**2 + (turbineY[j] - turbineY[i])**2 k += 1 - outputs['wtSeparationSquared'] = separation_squared -# print("wt sep") -# print(separation_squared) -# print() - + return separation_squared + + def compute_partials(self, inputs, J): - # obtain necessary inputs turbineX = inputs['turbineX'] turbineY = inputs['turbineY'] + dSdx, dSdy = self._compute_partials(turbineX, turbineY) + # populate Jacobian dict + J['wtSeparationSquared', 'turbineX'] = dSdx + J['wtSeparationSquared', 'turbineY'] = dSdy + + + def _compute_partials(self, turbineX, turbineY): # get number of turbines nTurbines = self.nTurbines # initialize gradient calculation array -# dS = np.zeros((int((nTurbines - 1.) * nTurbines / 2.), 2 * nTurbines)) dSdx = np.zeros((int((nTurbines - 1.) * nTurbines / 2.), nTurbines)) #col: dx_1-dx_n, row: d12, d13,..,d1n, d23..d2n,.. dSdy = np.zeros((int((nTurbines - 1.) * nTurbines / 2.), nTurbines)) @@ -71,20 +76,13 @@ class SpacingComp(ExplicitComponent): for i in range(0, nTurbines): for j in range(i + 1, nTurbines): # separation wrt Xj -# dS[k, j] = 2 * (turbineX[j] - turbineX[i]) dSdx[k,j]= 2 * (turbineX[j] - turbineX[i]) # separation wrt Xi -# dS[k, i] = -2 * (turbineX[j] - turbineX[i]) dSdx[k, i] = -2 * (turbineX[j] - turbineX[i]) # separation wrt Yj -# dS[k, j + nTurbines] = 2 * (turbineY[j] - turbineY[i]) dSdy[k, j] = 2 * (turbineY[j] - turbineY[i]) # separation wrt Yi -# dS[k, i + nTurbines] = -2 * (turbineY[j] - turbineY[i]) dSdy[k, i] = -2 * (turbineY[j] - turbineY[i]) # increment turbine pair counter k += 1 - - # populate Jacobian dict - J['wtSeparationSquared', 'turbineX'] = dSdx - J['wtSeparationSquared', 'turbineY'] = dSdy + return dSdx, dSdy diff --git a/topfarm/utils.py b/topfarm/utils.py index 16cdaaf76a3f39cc1eba6099b6be60400d5c87eb..f5256c7e1e4e62b6c0ffccfceab6859c34ba12a1 100644 --- a/topfarm/utils.py +++ b/topfarm/utils.py @@ -1,14 +1,18 @@ +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 -from topfarm.constraint_components.boundary_component import \ - PolygonBoundaryComp -from topfarm.constraint_components.spacing_component import SpacingComp def pos_from_case(case_recorder_filename): + ''' + Input: a recorded optimization case + Returns the positions at the last step of that optimization + ''' if not os.path.exists(case_recorder_filename): string = 'The specified optimization recording does not exists: ' string += case_recorder_filename @@ -25,6 +29,10 @@ def pos_from_case(case_recorder_filename): def latest_id(case_recorder_dir): + ''' + Input: path to the directory of recorded optimizations + Returns the absolute path to the most recent recording in that directory + ''' files = os.listdir(case_recorder_dir) files = [x for x in files if x.startswith('cases_') and x.endswith('.sql')] if len(files) == 0: @@ -40,32 +48,33 @@ def latest_id(case_recorder_dir): def random_positions(boundary, n_wt, n_iter, step_size, min_space, plot=False, verbose=True): - ba = np.array(boundary).T - x_max = max(ba[0]) - x_min = min(ba[0]) - y_max = max(ba[1]) - y_min = min(ba[1]) - x = np.random.rand(n_wt) - y = np.random.rand(n_wt) - xm = (x_max+x_min)/2 - ym = (y_max+y_min)/2 - xa = (x_max-xm) - ya = (y_max-ym) - x *= 2*xa - y *= 2*ya - x += xm-xa - y += ym-ya + ''' + Input: + boundary: list of tuples, e.g.: [(0, 0), (6, 1), (7, -11), (-1, -10)] + n_wt: number of wind turbines + n_iter: number of iterations allowed to try and satisfy the minimum + spacing constraint + step_size: the multiplier on the gradient that the turbines are moved + in each step + min_space: the minimum spacing between turbines + 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, plot, verbose) + return turbines + +def _random_positions(x, y, boundary, n_wt, n_iter, step_size, min_space, + plot, verbose): boundary_comp = PolygonBoundaryComp(boundary, n_wt) spacing_comp = SpacingComp(nTurbines=n_wt) turbineX = copy(x) turbineY = copy(y) - inputs = {} - outputs = {} - J = {} - inputs['turbineX'] = turbineX - inputs['turbineY'] = turbineY 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: @@ -73,11 +82,8 @@ def random_positions(boundary, n_wt, n_iter, step_size, min_space, plt.plot(bx, by) plt.plot(x, y, '.') for j in range(n_iter): - spacing_comp.compute(inputs, outputs) - spacing_comp.compute_partials(inputs, J) - dx = J['wtSeparationSquared', 'turbineX'] - dy = J['wtSeparationSquared', 'turbineY'] - dist = outputs['wtSeparationSquared'] + dist = spacing_comp._compute(turbineX, turbineY) + dx, dy = spacing_comp._compute_partials(turbineX, turbineY) index = np.argmin(dist) if dist[index] < min_space2: turbineX += dx[index]*step_size @@ -91,8 +97,7 @@ def random_positions(boundary, n_wt, n_iter, step_size, min_space, break if plot: plt.plot(turbineX, turbineY, 'o') - spacing_comp.compute(inputs, outputs) - dist = outputs['wtSeparationSquared'] + dist = spacing_comp._compute(turbineX, turbineY) index = np.argmin(dist) if verbose: print('Spacing obtained: {}'.format(dist[index]**0.5)) @@ -101,6 +106,10 @@ def random_positions(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): for i in range(0, n_wt): dng = boundary_comp.calc_distance_and_gradients(turbineX, @@ -111,21 +120,19 @@ def _contain(n_wt, turbineX, turbineY, boundary_comp): dy = dng[2][i] turbineX[i] -= dx*dist*1.1 turbineY[i] -= dy*dist*1.1 - dng = boundary_comp.calc_distance_and_gradients(turbineX, - turbineY) - dist = dng[0][i] return turbineX, turbineY if __name__ == '__main__': - crf = r"C:\Sandbox\Git\TopFarm2\topfarm\cases_20180621_111710.sql" - case_recorder_filename = crf - turbines = pos_from_case(case_recorder_filename) - print(turbines) - - case_recorder_dir = r'C:\Sandbox\Git\TopFarm2\topfarm' - latest_id = latest_id(case_recorder_dir) - print(latest_id) +# this_dir = os.path.dirname(os.path.abspath(__file__)) +# crf = r"C:\Sandbox\Git\TopFarm2\topfarm\cases_20180621_111710.sql" +# case_recorder_filename = crf +# turbines = pos_from_case(case_recorder_filename) +# print(turbines) +# +# case_recorder_dir = r'C:\Sandbox\Git\TopFarm2\topfarm' +# latest_id = latest_id(case_recorder_dir) +# print(latest_id) boundary = [(0, 0), (6, 1), (7, -11), (-1, -10)] n_wt = 20