diff --git a/topfarm/_topfarm.py b/topfarm/_topfarm.py
index d3de26ec9fef144aad7691a50fa8481faaec62e3..ac9257d55a6d759ac2c8db4d117190bb06174083 100644
--- a/topfarm/_topfarm.py
+++ b/topfarm/_topfarm.py
@@ -1,6 +1,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 smart_start
 import os
 import time
 import numpy as np
@@ -240,6 +241,11 @@ class TurbineXYZOptimizationProblem(TopFarmProblem):
     def z_boundary(self):
         return self.boundary_comp.z_boundary
 
+    def smart_start(self, x, y ,val):
+        x, y = smart_start(x, y ,val, self.n_wt, self.min_spacing)
+        self.update_state({'turbineX': x, 'turbineY': y})
+        return x, y
+
 
 class InitialXYZOptimizationProblem(TurbineXYZOptimizationProblem):
     def __init__(self, cost_comp, turbineXYZ, boundary_comp=None, min_spacing=None,
diff --git a/topfarm/tests/test_topfarm_utils/test_smart_start.py b/topfarm/tests/test_topfarm_utils/test_smart_start.py
new file mode 100644
index 0000000000000000000000000000000000000000..eeae4dac18206e8c051df93d3de74efdd60776c9
--- /dev/null
+++ b/topfarm/tests/test_topfarm_utils/test_smart_start.py
@@ -0,0 +1,60 @@
+'''
+Created on 17. jul. 2018
+
+@author: mmpe
+'''
+from topfarm.utils import smart_start
+import numpy as np
+
+x = np.arange(0, 20, 0.1)
+y = np.arange(0, 10, 0.1)
+YY, XX = np.meshgrid(y, x)
+val = np.sin(XX) + np.sin(YY)
+N_WT = 20
+min_space = 2.1
+xs_ref = [1.6,
+ 14.100000000000001,
+ 1.6,
+ 7.9,
+ 14.100000000000001,
+ 7.9,
+ 19.900000000000002,
+ 19.900000000000002,
+ 5.800000000000001,
+ 7.800000000000001,
+ 14.200000000000001,
+ 1.5,
+ 5.800000000000001,
+ 16.2,
+ 16.2,
+ 1.6,
+ 3.7,
+ 14.100000000000001,
+ 3.7,
+ 7.9]
+ys_ref = [1.6,
+ 1.6,
+ 7.9,
+ 1.6,
+ 7.9,
+ 7.9,
+ 1.6,
+ 7.9,
+ 7.800000000000001,
+ 5.800000000000001,
+ 5.800000000000001,
+ 5.800000000000001,
+ 1.5,
+ 7.800000000000001,
+ 1.5,
+ 3.7,
+ 1.6,
+ 3.7,
+ 7.9,
+ 3.7]
+
+
+def tests_smart_start():
+    xs, ys = smart_start(XX, YY, val, N_WT, min_space)
+    return np.testing.assert_allclose(np.array([xs, ys]),
+                                      np.array([xs_ref, ys_ref]))
diff --git a/topfarm/utils.py b/topfarm/utils.py
index a6847a55492fc264f29fe056e82c046efa723a96..f798844258663c15707fa63293f79c16edc6063a 100644
--- a/topfarm/utils.py
+++ b/topfarm/utils.py
@@ -4,10 +4,6 @@ from topfarm.constraint_components.spacing_component import SpacingComp
 import os
 import numpy as np
 from openmdao.api import CaseReader
-import matplotlib.pyplot as plt
-
-
-
 
 
 def pos_from_case(case_recorder_filename):
@@ -89,6 +85,7 @@ def _shuffle_postions_rel(init_pos, offset, boundary, n_wt, plot):
     ba = np.array(boundary).T
     turbines = np.array(init_pos) + np.random.rand(n_wt, 2)*2*offset-offset
     if plot:
+        plt.figure()
         plt.cla()
         plt.plot(ba[0], ba[1])
         plt.plot(turbineX, turbineY, '.')
@@ -103,6 +100,7 @@ def _shuffle_positions_abs(turbineX, turbineY, boundary, n_wt, n_iter,
     min_space2 = min_space**2
     ba = np.array(boundary).T
     if plot:
+        plt.figure()
         plt.cla()
         plt.plot(ba[0], ba[1])
         plt.plot(turbineX, turbineY, '.')
@@ -144,18 +142,42 @@ def _move_inside_boundary(n_wt, turbineX, turbineY, boundary_comp, pad):
     return turbineX, turbineY
 
 
-if __name__ == '__main__':
-    this_dir = os.getcwd()
-    crf = r"tests\test_files\recordings\cases_20180703_152607.sql"
-    case_recorder_filename = crf
-    path = os.path.join(this_dir, crf)
-    turbines = pos_from_case(path)
-    print(turbines)
+def smart_start(x, y, val, N_WT, min_space):
+    '''
+    Selects the a number of gridpoints (N_WT) in the grid defined by x and y,
+    where val has the maximum value, while chosen points spacing (min_space)
+    is respected.
+    Input:
+        x and y: arrays (nD) with coordinates
+        val: array (nD), the corresponding value of the desired variable in the
+        grid points. This could be e.g. the AEP or wind speed.
+        N_WT: integer, number of wind turbines
+        min_space: float, minimum space between turbines
+    Output:
+        Positions where the aep or wsp is highest, while
+        respecting the minimum spacing requirement.
+    '''
+    arr = np.array([x.flatten(), y.flatten(), val.flatten()])
+    xs, ys = [], []
+    for i in range(N_WT):
+        try:
+            max_ind = np.argmax(arr[2])
+            x0 = arr[0][max_ind]
+            y0 = arr[1][max_ind]
+            xs.append(x0)
+            ys.append(y0)
+            index = np.where((arr[0]-x0)**2+(arr[1]-y0)**2 >= min_space**2)[0]
+            arr = arr[:, index]
+        except ValueError:
+            xs.append(np.nan)
+            ys.append(np.nan)
+            print('Could not respect the spacing constraint')
+    return xs, ys
 
-    case_recorder_dir = r'tests\test_files\recordings'
-    latest_id = latest_id(case_recorder_dir)
-    print(latest_id)
 
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    plt.clf()
     boundary = [(0, 0), (6, 1), (7, -11), (-1, -10), (0, 0)]
     n_wt = 20
     init_pos = np.column_stack((np.random.randint(0, 6, (n_wt)),
@@ -164,3 +186,22 @@ if __name__ == '__main__':
     turbines = shuffle_positions(boundary, n_wt, min_space, init_pos,
                                  shuffle_type='rel')
     print(turbines)
+
+    x = np.arange(0, 20, 0.1)
+    y = np.arange(0, 10, 0.1)
+    YY, XX = np.meshgrid(y, x)
+    val = np.sin(XX) + np.sin(YY)
+    N_WT = 30
+    min_space = 2.1
+    xs, ys = smart_start(XX, YY, val, N_WT, min_space)
+    plt.figure(1)
+    c = plt.contourf(XX, YY, val, 100)
+    plt.colorbar(c)
+    for i in range(N_WT):
+        circle = plt.Circle((xs[i], ys[i]), min_space/2, color='b', fill=False)
+        plt.gcf().gca().add_artist(circle)
+        plt.plot(xs[i], ys[i], 'rx')
+    plt.axis('equal')
+    plt.show()
+
+