diff --git a/topfarm/_topfarm.py b/topfarm/_topfarm.py
index d3de26ec9fef144aad7691a50fa8481faaec62e3..75d62a657b96b5f710787513a3c0d136b804ff86 100644
--- a/topfarm/_topfarm.py
+++ b/topfarm/_topfarm.py
@@ -15,17 +15,31 @@ from openmdao.api import Problem, ScipyOptimizeDriver, IndepVarComp
class TopFarmProblem(Problem):
- def __init__(self, cost_comp, record_id, expected_cost):
+ def __init__(self, cost_comp, driver, plot_comp, record_id, expected_cost):
Problem.__init__(self)
if isinstance(cost_comp, TopFarmProblem):
cost_comp = cost_comp.as_component()
cost_comp.parent = self
- self.record_id = record_id
self.cost_comp = cost_comp
+
+ if isinstance(driver, list):
+ driver = DOEDriver(ListGenerator(driver))
+ elif isinstance(driver, DOEGenerator):
+ driver = DOEDriver(generator=driver)
+ self.driver = driver
+
+ self.plot_comp = plot_comp
+
+ self.record_id = record_id
+
self.indeps = self.model.add_subsystem('indeps', IndepVarComp(), promotes=['*'])
self.model.add_subsystem('cost_comp', cost_comp, promotes=['*'])
self.model.add_objective('cost', scaler=1 / abs(expected_cost))
+ if plot_comp:
+ self.model.add_subsystem('plot_comp', plot_comp, promotes=['*'])
+ plot_comp.problem = self
+
@property
def cost(self):
return self['cost'][0]
@@ -54,10 +68,10 @@ class TopFarmProblem(Problem):
def evaluate(self, state={}):
t = time.time()
self.update_state(state)
- rec = ListRecorder()
- self.driver.add_recorder(rec)
+ self.recorder = ListRecorder()
+ self.driver.add_recorder(self.recorder)
self.run_model()
- self.driver._rec_mgr._recorders.remove(rec)
+ self.driver._rec_mgr._recorders.remove(self.recorder)
print("Evaluated in\t%.3fs" % (time.time() - t))
return self.cost, self.state
@@ -80,8 +94,8 @@ class TopFarmProblem(Problem):
try:
self.update_state({k: self.recorder.get(k)[-1] for k in self.state.keys() if k not in state})
except ValueError:
- pass # loaded state does not fit into dimension of current state
-
+ pass # loaded state does not fit into dimension of current state
+
self.driver.add_recorder(self.recorder)
self.run_driver()
self.cleanup()
@@ -126,8 +140,8 @@ class TopFarmProblem(Problem):
class TurbineTypeOptimizationProblem(TopFarmProblem):
- def __init__(self, cost_comp, turbineTypes, lower, upper, driver, record_id=None, expected_cost=1):
- TopFarmProblem.__init__(self, cost_comp, record_id, expected_cost)
+ def __init__(self, cost_comp, turbineTypes, lower, upper, driver, plot_comp=None, record_id=None, expected_cost=1):
+ TopFarmProblem.__init__(self, cost_comp, driver, plot_comp, record_id, expected_cost)
self.turbineTypes = turbineTypes
self.lower = lower
self.upper = upper
@@ -141,13 +155,15 @@ class TurbineTypeOptimizationProblem(TopFarmProblem):
self.model.add_constraint('turbineType', lower=lim[:, 0], upper=lim[:, 1])
self.indeps.add_output('turbineType', np.array(turbineTypes).astype(np.int))
- self.driver = driver
+
self.model.add_design_var('turbineType', lower=lim[:, 0], upper=lim[:, 1])
+ if self.plot_comp:
+ plot_comp.n_wt = n_wt
self.setup(check=True, mode='fwd')
@property
def state(self):
- state = {'turbineType': np.array(self['turbineType']).astype(np.int)}
+ state = {'turbineType': np.round(self['turbineType']).astype(np.int)}
state.update(TopFarmProblem.state.fget(self))
return state
@@ -155,11 +171,16 @@ class TurbineTypeOptimizationProblem(TopFarmProblem):
class TurbineXYZOptimizationProblem(TopFarmProblem):
def __init__(self, cost_comp, turbineXYZ, boundary_comp, min_spacing=None,
driver=ScipyOptimizeDriver(), plot_comp=None, record_id=None, expected_cost=1):
- self.turbineXYZ = turbineXYZ
- TopFarmProblem.__init__(self, cost_comp, record_id, expected_cost)
-
turbineXYZ = np.array(turbineXYZ)
+ self.turbineXYZ = turbineXYZ
self.n_wt = n_wt = turbineXYZ.shape[0]
+
+ if plot_comp:
+ if plot_comp == "default":
+ plot_comp = PlotComp()
+
+ TopFarmProblem.__init__(self, cost_comp, driver, plot_comp, record_id, expected_cost)
+
turbineXYZ = np.hstack((turbineXYZ, np.zeros((n_wt, 4 - turbineXYZ.shape[1]))))
self.min_spacing = min_spacing
@@ -169,7 +190,7 @@ class TurbineXYZOptimizationProblem(TopFarmProblem):
self.boundary_comp = boundary_comp
boundary_comp.setup_as_constraints(self)
- do = driver.options
+ do = self.driver.options
if len(boundary_comp.xy_boundary) > 0:
ref0_x, ref0_y = self.boundary_comp.xy_boundary.min(0)
@@ -202,20 +223,13 @@ class TurbineXYZOptimizationProblem(TopFarmProblem):
self.indeps.add_output('turbineY', turbineXYZ[:, 1], units='m')
self.indeps.add_output('turbineZ', turbineXYZ[:, 2], units='m')
- self.driver = driver
-
if plot_comp:
- if plot_comp == "default":
- plot_comp = PlotComp()
-
plot_comp.n_wt = n_wt
if self.boundary_comp:
plot_comp.n_vertices = len(self.boundary_comp.xy_boundary)
else:
plot_comp.n_vertices = 0
- self.model.add_subsystem('plot_comp', plot_comp, promotes=['*'])
- self.plot_comp = plot_comp
self.setup(check=True, mode='fwd')
@property
@@ -246,10 +260,6 @@ class InitialXYZOptimizationProblem(TurbineXYZOptimizationProblem):
driver=None, plot_comp=None):
# if driver is None:
# driver = DOEDriver(shuffle_generator(self, 10))
- if isinstance(driver, list):
- driver = DOEDriver(ListGenerator(driver))
- elif isinstance(driver, DOEGenerator):
- driver = DOEDriver(generator=driver)
TurbineXYZOptimizationProblem.__init__(self, cost_comp, turbineXYZ,
boundary_comp, min_spacing,
driver=driver, plot_comp=plot_comp)
diff --git a/topfarm/constraint_components/constrained_generator.py b/topfarm/constraint_components/constrained_generator.py
index 612944491e8c3901bdfe73cb931588b7bc06785d..cecaaaf1f84c57d5446943683a6c39837a787342 100644
--- a/topfarm/constraint_components/constrained_generator.py
+++ b/topfarm/constraint_components/constrained_generator.py
@@ -2,6 +2,30 @@ from openmdao.drivers.doe_generators import DOEGenerator, ListGenerator
import numpy as np
+
+class ConstrainedDiscardXYZGenerator(DOEGenerator):
+ def __init__(self, generator):
+ DOEGenerator.__init__(self)
+ if isinstance(generator, list):
+ generator = ListGenerator(generator)
+ self.generator = generator
+
+ def __call__(self, design_vars, model=None):
+ xy_boundary_comp = model._get_subsystem('xy_bound_comp')
+ spacing_comp = model._get_subsystem('spacing_comp')
+ for xyz_tuple_lst in self.generator(design_vars, model=model):
+ x, y, z = ([np.array(t[1]) for t in xyz_tuple_lst] + [0])[:3]
+ if spacing_comp:
+ dist = spacing_comp._compute(x, y)
+ if np.any(dist < spacing_comp.min_spacing**2):
+ continue
+
+ if xy_boundary_comp:
+ dist_to_boundary =xy_boundary_comp.distances(x,y)
+ if np.any(dist_to_boundary<0):
+ continue
+ yield xyz_tuple_lst
+
class ConstrainedXYZGenerator(DOEGenerator):
def __init__(self, generator, n_iter=1000, step_size=0.1, offset=0.5, verbose=False):
DOEGenerator.__init__(self)
@@ -17,16 +41,15 @@ class ConstrainedXYZGenerator(DOEGenerator):
xy_boundary_comp = model._get_subsystem('xy_bound_comp')
spacing_comp = model._get_subsystem('spacing_comp')
for xyz_tuple_lst in self.generator(design_vars, model=model):
- x, y, z = [np.array(t[1]) for t in xyz_tuple_lst]
+ x, y, z = ([np.array(t[1]).astype(np.float) for t in xyz_tuple_lst] + [0])[:3]
if spacing_comp:
for j in range(self.n_iter):
dist = spacing_comp._compute(x, y)
dx, dy = spacing_comp._compute_partials(x, y)
- dx, dy = dx[0], dy[0]
- index = np.argmin(dist)
+ index = int(np.argmin(dist))
if dist[index] < spacing_comp.min_spacing**2 or j == 0:
- x[index] += dx[index] * self.step_size
- y[index] += dy[index] * self.step_size
+ x += dx[index] * self.step_size
+ y += dy[index] * self.step_size
if xy_boundary_comp:
x, y, z = xy_boundary_comp.move_inside(x, y, z)
else:
diff --git a/topfarm/plotting.py b/topfarm/plotting.py
index 0f710f2f5977d5aaf936a5d7279005e4cea1c1be..4474ccad01e2d63cdd24df30fbbac1935271cf3f 100644
--- a/topfarm/plotting.py
+++ b/topfarm/plotting.py
@@ -21,15 +21,19 @@ def mypause(interval):
class PlotComp(ExplicitComponent):
colors = ['b', 'r', 'm', 'c', 'g', 'y', 'orange', 'indigo', 'grey'] * 100
- def __init__(self, memory=10, delay=0.001, plot_initial=True,
- animate=False):
+ def __init__(self, memory=10, delay=0.001, plot_initial=True, ax=None):
ExplicitComponent.__init__(self)
+ self.ax_ = ax
self.memory = memory
self.delay = delay
self.plot_initial = plot_initial
self.history = []
self.counter = 0
- self.animate = animate
+ self.by_pass = False
+
+ @property
+ def ax(self):
+ return self.ax_ or plt.gca()
def show(self):
plt.show()
@@ -38,66 +42,92 @@ class PlotComp(ExplicitComponent):
self.add_input('turbineX', np.zeros(self.n_wt), units='m')
self.add_input('turbineY', np.zeros(self.n_wt), units='m')
self.add_input('cost', 0.)
- self.add_input('boundary', np.zeros((self.n_vertices, 2)), units='m')
-
- def init_plot(self, boundary):
- plt.cla()
- plt.axis('equal')
- b = np.r_[boundary[:], boundary[:1]]
- plt.plot(b[:, 0], b[:, 1], 'k')
- mi = b.min(0)
- ma = b.max(0)
+ if hasattr(self, 'n_vertices'):
+ self.add_input('boundary', np.zeros((self.n_vertices, 2)), units='m')
+
+ def init_plot(self, limits):
+ self.ax.cla()
+ self.ax.axis('equal')
+ mi = limits.min(0)
+ ma = limits.max(0)
ra = ma - mi
ext = .1
xlim, ylim = np.array([mi - ext * ra, ma + ext * ra]).T
- plt.xlim(xlim)
- plt.ylim(ylim)
+ self.ax.set_xlim(xlim)
+ self.ax.set_ylim(ylim)
+
+ def plot_boundary(self, boundary):
+ b = np.r_[boundary[:], boundary[:1]]
+ plt.plot(b[:, 0], b[:, 1], 'k')
+
+ def plot_history(self, x, y):
+ rec = self.problem.recorder
+ if rec.num_cases > 0:
+ x = np.r_[rec['turbineX'][-self.memory:], [x]]
+ y = np.r_[rec['turbineY'][-self.memory:], [y]]
+ for c, x_, y_ in zip(self.colors, x.T, y.T):
+ self.ax.plot(x_, y_, '--', color=c)
+
+ def plot_initial2current(self, x0, y0, x, y):
+ rec = self.problem.recorder
+ if rec.num_cases > 0:
+ x0 = rec['turbineX'][0]
+ y0 = rec['turbineY'][0]
+ for c, x0_, y0_, x_, y_ in zip(self.colors, x0, y0, x, y):
+ self.ax.plot((x0_, x_), (y0_, y_), '-', color=c)
+
+ def plot_current_position(self, x, y):
+ for c, x_, y_ in zip(self.colors, x, y):
+ self.ax.plot(x_, y_, 'o', color=c, ms=5)
+ self.ax.plot(x_, y_, 'xk', ms=4)
+
+ def set_title(self, cost0, cost):
+ rec = self.problem.recorder
+ if cost0 != 0:
+ self.ax.set_title("%d: %f (%.2f%%)" % (rec.num_cases, cost, (cost0 - cost) / cost0 * 100))
+ else:
+ self.ax.set_title("%d: %f" % (rec.num_cases, cost))
+
+ def get_initial(self):
+ rec = self.problem.recorder
+ if rec.num_cases > 0:
+ x0 = rec['turbineX'][0]
+ y0 = rec['turbineY'][0]
+ cost0 = rec['cost'][0]
+ return x0, y0, cost0
def compute(self, inputs, outputs):
- x = inputs['turbineX']
- y = inputs['turbineY']
- cost = inputs['cost'][0]
- if not hasattr(self, "initial"):
- self.initial = np.array([x, y]).T, cost
- self.history = [(x.copy(), y.copy())] + self.history[:self.memory]
-
- boundary = inputs['boundary']
- self.init_plot(boundary)
- plt.title("%d: %f (%.2f%%)" % (self.counter, cost,
- (self.initial[1]-cost)/self.initial[1]*100))
-
- history_arr = np.array(self.history)
- for i, c, x_, y_ in zip(range(len(x)), self.colors, x, y):
- if self.plot_initial:
- plt.plot([self.initial[0][i, 0], x_],
- [self.initial[0][i, 1], y_], '-', color=c, lw=1)
- plt.plot(history_arr[:, 0, i], history_arr[:, 1, i], '.--',
- color=c, lw=1)
- plt.plot(x_, y_, 'o', color=c, ms=5)
- plt.plot(x_, y_, 'x' + 'k', ms=4)
- self.temp = '../__animations__'
- if not os.path.exists(self.temp):
- os.makedirs(self.temp)
- if self.animate:
- path = os.path.join(self.temp,
- 'plot_{:05d}.png'.format(self.counter))
- plt.savefig(path)
-
- if self.counter == 0:
- plt.pause(.01)
- mypause(self.delay)
-
- self.counter += 1
-
- def run_animate(self, anim_time=10, verbose=False):
- N = anim_time/self.counter
- string = 'ffmpeg -f image2 -r 1/'
- string += '{} -i {}//plot_%05d.png'.format(N, self.temp)
- string += ' -vcodec mpeg4 -y {}//animation.mp4'.format(self.temp)
- if verbose:
- print('\nCreating animation:')
- print(string)
- os.system(string)
+ if self.by_pass is False:
+ x = inputs['turbineX']
+ y = inputs['turbineY']
+ cost = inputs['cost'][0]
+
+ if 'boundary' in inputs:
+ boundary = inputs['boundary']
+ self.init_plot(boundary)
+ self.plot_boundary(boundary)
+ else:
+ self.init_plot(np.array([x, y]).T)
+
+ initial = self.get_initial()
+
+ if initial is not None:
+ x0, y0, cost0 = initial
+ if self.plot_initial:
+ self.plot_initial2current(x0, y0, x, y)
+ if self.memory > 0:
+ self.plot_history(x, y)
+ else:
+ cost0 = cost
+
+ self.plot_current_position(x, y)
+ self.set_title(cost0, cost)
+
+ if self.counter == 0:
+ plt.pause(.01)
+ mypause(self.delay)
+
+ self.counter += 1
class NoPlot(PlotComp):
@@ -108,10 +138,35 @@ class NoPlot(PlotComp):
pass
def setup(self):
- self.add_input('turbineX', np.zeros(self.n_wt), units='m')
- self.add_input('turbineY', np.zeros(self.n_wt), units='m')
self.add_input('cost', 0.)
- self.add_input('boundary', np.zeros((self.n_vertices, 2)), units='m')
def compute(self, inputs, outputs):
pass
+
+
+class TurbineTypePlotComponent(PlotComp):
+ colors = ['b', 'r', 'm', 'c', 'g', 'y', 'orange', 'indigo', 'grey'] * 100
+ markers = np.array(list(".ov^<>12348spP*hH+xXDd|_"))
+
+ def __init__(self, turbine_type_names, memory=10, delay=0.001, plot_initial=True):
+ self.turbine_type_names = turbine_type_names
+ PlotComp.__init__(self, memory=memory, delay=delay, plot_initial=plot_initial)
+
+ def setup(self):
+ PlotComp.setup(self)
+ self.add_input('turbineType', np.zeros(self.n_wt, dtype=np.int))
+
+ def compute(self, inputs, outputs):
+ self.types = np.asarray(inputs['turbineType'], dtype=np.int)
+ PlotComp.compute(self, inputs, outputs)
+
+ def init_plot(self, limits):
+ PlotComp.init_plot(self, limits)
+ for m, n in zip(self.markers, self.turbine_type_names):
+ self.ax.plot([], [], m + 'k', label=n)
+ self.ax.legend()
+
+ def plot_current_position(self, x, y):
+ for m, c, x_, y_ in zip(self.markers[self.types], self.colors, x, y):
+ #self.ax.plot(x_, y_, 'o', color=c, ms=5)
+ self.ax.plot(x_, y_, m + 'k', markeredgecolor=c, markeredgewidth=1, ms=8)
diff --git a/topfarm/recorders.py b/topfarm/recorders.py
index 216dcc1ccb80f9aac158a2a175dd2bdbbf3c13dc..03252a43afdc2774925d3ff29859a689f4c46183 100644
--- a/topfarm/recorders.py
+++ b/topfarm/recorders.py
@@ -84,6 +84,10 @@ class ListRecorder(BaseRecorder):
self._abs2meta[name]['explicit'] = True
if name in states:
self._abs2meta[name]['explicit'] = False
+
+ @property
+ def num_cases(self):
+ return len(self.driver_iteration_lst)
@property
def driver_cases(self):
diff --git a/topfarm/tests/test_topfarm_utils/test_parellel_runner.py b/topfarm/tests/test_topfarm_utils/test_parellel_runner.py
index f14439207c87c6f3a2d30172515033e8695be7d4..4a44a4c94efa0f9b5108a17a618db9774ffe2911 100644
--- a/topfarm/tests/test_topfarm_utils/test_parellel_runner.py
+++ b/topfarm/tests/test_topfarm_utils/test_parellel_runner.py
@@ -1,12 +1,10 @@
from openmdao.drivers.doe_driver import DOEDriver
-from openmdao.drivers.doe_generators import FullFactorialGenerator,\
- ListGenerator, UniformGenerator
+from openmdao.drivers.doe_generators import UniformGenerator
import pytest
from topfarm._topfarm import InitialXYZOptimizationProblem
import numpy as np
from topfarm.cost_models.dummy import DummyCost
-from topfarm.tests import npt, uta
-from topfarm.constraint_components.constrained_generator import ConstrainedXYZGenerator
+from topfarm.tests import npt
from topfarm.constraint_components.boundary_component import BoundaryComp
from topfarm.parallel_runner import ParallelRunner
diff --git a/topfarm/tests/test_with_dummy.py b/topfarm/tests/test_with_dummy.py
index d9563e79bc193612c8a717d640540cdbd4cb76b9..103247584beb7845249c9db8dc8a7a9b68978cb0 100644
--- a/topfarm/tests/test_with_dummy.py
+++ b/topfarm/tests/test_with_dummy.py
@@ -6,52 +6,49 @@ import unittest
import pytest
import numpy as np
from topfarm.cost_models.dummy import DummyCost, DummyCostPlotComp
-
-
-class Test(unittest.TestCase): # unittest version
-
- def test_optimize_4tb(self):
- """Optimize 4-turbine layout and check final positions
-
- The farm boundaries and min spacing are chosen such that the desired
- turbine positions are not within the boundaries or constraints.
- """
-
- # test options
- dec_prec = 4 # decimal precision for comparison
-
- # given
- boundary = [(0, 0), (6, 0), (6, -10), (0, -10)] # turbine boundaries
- initial = [[6, 0], [6, -8], [1, 1], [-1, -8]] # initial turbine layouts
- desired = [[3, -3], [7, -7], [4, -3], [3, -7]] # desired turbine layouts
- optimal = np.array([[2.5, -3], [6, -7], [4.5, -3], [3, -7]]) # optimal turbine layout
- min_spacing = 2 # min distance between turbines
-
- # when
- tf = TopFarm(initial, DummyCost(desired, ['turbineX','turbineY']), min_spacing,
- boundary=boundary, record_id=None)
- tf.optimize()
- tb_pos = tf.turbine_positions
-
- # then
- tol = 1e-6
- self.assertGreater(sum((tb_pos[2] - tb_pos[0])**2), 2**2 - tol) # check min spacing
- self.assertLess(tb_pos[1][0], 6 + tol) # check within border
- np.testing.assert_array_almost_equal(tb_pos[:,:2], optimal, dec_prec)
-
- def testDummyCostPlotComp(self):
- if os.name == 'posix' and "DISPLAY" not in os.environ:
- pytest.xfail("No display")
-
- desired = [[3, -3], [7, -7], [4, -3], [3, -7]]
- tf = TopFarm(turbines=[[6, 0], [6, -8], [1, 1], [-1, -8]],
- cost_comp=DummyCost(desired, ['turbineX','turbineY']),
- min_spacing=2,
- boundary=[(0, 0), (6, 0), (6, -10), (0, -10)],
- plot_comp = DummyCostPlotComp(desired),
- record_id=None)
- tf.evaluate()
-
-
-if __name__ == "__main__":
- unittest.main()
+from topfarm.tests import uta
+
+
+def test_optimize_4tb():
+ """Optimize 4-turbine layout and check final positions
+
+ The farm boundaries and min spacing are chosen such that the desired
+ turbine positions are not within the boundaries or constraints.
+ """
+
+ # test options
+ dec_prec = 4 # decimal precision for comparison
+
+ # given
+ boundary = [(0, 0), (6, 0), (6, -10), (0, -10)] # turbine boundaries
+ initial = [[6, 0], [6, -8], [1, 1], [-1, -8]] # initial turbine layouts
+ desired = [[3, -3], [7, -7], [4, -3], [3, -7]] # desired turbine layouts
+ optimal = np.array([[2.5, -3], [6, -7], [4.5, -3], [3, -7]]) # optimal turbine layout
+ min_spacing = 2 # min distance between turbines
+
+ # when
+ tf = TopFarm(initial, DummyCost(desired, ['turbineX', 'turbineY']), min_spacing,
+ boundary=boundary, record_id=None)
+ tf.optimize()
+ tb_pos = tf.turbine_positions
+
+ # then
+ tol = 1e-6
+ uta.assertGreater(sum((tb_pos[2] - tb_pos[0])**2), 2**2 - tol) # check min spacing
+ uta.assertLess(tb_pos[1][0], 6 + tol) # check within border
+ np.testing.assert_array_almost_equal(tb_pos[:, :2], optimal, dec_prec)
+
+
+def testDummyCostPlotComp():
+ if os.name == 'posix' and "DISPLAY" not in os.environ:
+ pytest.xfail("No display")
+
+ desired = [[3, -3], [7, -7], [4, -3], [3, -7]]
+ tf = TopFarm(turbines=[[6, 0], [6, -8], [1, 1], [-1, -8]],
+ cost_comp=DummyCost(desired, ['turbineX', 'turbineY']),
+ min_spacing=2,
+ boundary=[(0, 0), (6, 0), (6, -10), (0, -10)],
+ plot_comp=DummyCostPlotComp(desired),
+ record_id=None)
+ tf.evaluate()
+ tf.optimize()