updated utils script for clarity and simplicity and easier testing

updated spacing_component to avoid having to use openmdao compliant call to calculate distances and gradients.

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

topfarm/utils.py

+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