Commit 44f49c3a authored by Mikkel Friis-Møller's avatar Mikkel Friis-Møller
Browse files

Om316

parent 4777ee24
Pipeline #30170 passed with stages
in 15 minutes and 34 seconds
......@@ -33,7 +33,7 @@ test_topfarm_windows: # name the job what we like
script: # runs on windows machine due to tag below
- conda init powershell
- "if (test-path $PROFILE.CurrentUserAllHosts) { & $PROFILE.CurrentUserAllHosts}"
- conda activate py36_openmdao26
- conda activate om3
- pip install --upgrade git+https://gitlab.windenergy.dtu.dk/TOPFARM/PyWake.git
- pip install -e .
- python -m pytest --cov-report term-missing:skip-covered --cov=topfarm --cov-config .coveragerc --ignore=topfarm/cost_models/fuga/Colonel
......
import numpy as np
from openmdao.api import view_model
from openmdao.api import n2
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site
from py_wake.aep_calculator import AEPCalculator
from topfarm.cost_models.economic_models.turbine_cost import economic_evaluation
......@@ -42,18 +42,18 @@ def main():
input_keys=['x', 'y'],
n_wt=n_wt,
cost_function=lambda x, y, **_: windFarmModel(x=x, y=y).aep().sum(['wd', 'ws']) * 10**6,
output_key="aep",
output_keys="aep",
output_unit="GWh",
objective=False,
output_val=np.zeros(n_wt))
output_vals=np.zeros(n_wt))
irr_comp = CostModelComponent(
input_keys=['aep'],
n_wt=n_wt,
cost_function=irr_func,
output_key="irr",
output_keys="irr",
output_unit="%",
objective=True,
income_model=True)
maximize=True)
group = TopFarmGroup([aep_comp, irr_comp])
problem = TopFarmProblem(
design_vars=dict(zip('xy', site.initial_position.T)),
......@@ -64,7 +64,7 @@ def main():
plot_comp=plot_comp)
cost, state, recorder = problem.optimize()
# problem.evaluate()
# view_model(problem, outfile='example_4_integrated_optimization_aep_and_irr_n2.html', show_browser=False)
n2(problem, outfile='example_4_integrated_optimization_aep_and_irr_n2.html', show_browser=False)
main()
import numpy as np
from openmdao.api import view_model
from openmdao.api import n2
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site
from py_wake.aep_calculator import AEPCalculator
from topfarm.cost_models.economic_models.dtu_wind_cm_main import economic_evaluation
......@@ -54,18 +54,18 @@ def main():
input_keys=['x', 'y'],
n_wt=n_wt,
cost_function=lambda x, y, **_: windFarmModel(x=x, y=y).aep().sum(['wd', 'ws']) * 10**6,
output_key="aep",
output_keys="aep",
output_unit="kWh",
objective=False,
output_val=np.zeros(n_wt))
output_vals=np.zeros(n_wt))
irr_comp = CostModelComponent(
input_keys=['aep'],
n_wt=n_wt,
cost_function=irr_func,
output_key="irr",
output_keys="irr",
output_unit="%",
objective=True,
income_model=True)
maximize=True)
group = TopFarmGroup([aep_comp, irr_comp])
problem = TopFarmProblem(
design_vars=dict(zip('xy', site.initial_position.T)),
......@@ -76,7 +76,7 @@ def main():
plot_comp=plot_comp)
cost, state, recorder = problem.optimize()
# problem.evaluate()
# view_model(problem, outfile='example_4_integrated_optimization_aep_and_irr_n2.html', show_browser=False)
# n2(problem, outfile='example_4_integrated_optimization_aep_and_irr_n2.html', show_browser=False)
main()
......@@ -3,7 +3,7 @@ from py_wake.deficit_models.noj import NOJ
from py_wake.aep_calculator import AEPCalculator
import matplotlib.pyplot as plt
import numpy as np
from openmdao.api import view_model, IndepVarComp
from openmdao.api import n2, IndepVarComp
from topfarm import TopFarmGroup, TopFarmProblem
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from topfarm.cost_models.cost_model_wrappers import CostModelComponent, AEPCostModelComponent
......@@ -70,7 +70,7 @@ def main():
def cost_func(AEP, **kwargs):
return AEP
cost_comp = CostModelComponent(input_keys=[('AEP', [0])], n_wt=n_wt, cost_function=cost_func,
output_key="aep", output_unit="kWh", objective=True, income_model=True, input_units=['kW*h'])
output_keys="aep", output_unit="kWh", objective=True, maximize=True, input_units=['kW*h'])
group = TopFarmGroup([aep_comp, cost_comp])
boundary = np.array([(900, 1000), (2300, 1000), (2300, 2100), (900, 2100)]) # turbine boundaries
prob = TopFarmProblem(
......@@ -128,7 +128,7 @@ def main():
differentiable = True
prob = setup_prob(differentiable)
# view_model(prob)
# n2(prob)
cost_init, state_init = prob.evaluate()
tic = time.time()
cost, state, recorder = prob.optimize()
......
......@@ -7,7 +7,7 @@ from topfarm.cost_models.dummy import DummyCost, DummyCostPlotComp
from topfarm.easy_drivers import EasyScipyOptimizeDriver
from topfarm.plotting import NoPlot
import time
from openmdao.api import view_model
from openmdao.api import n2
......@@ -34,7 +34,7 @@ def main():
comps = [CostModelComponent('xy', 4,
cost_function=lambda x, y, i=i:wt_cost(i, x, y),
objective=False,
output_key='cost%d' % i) for i in range(n_wt)]
output_keys='cost%d' % i) for i in range(n_wt)]
def sum_map(**kwargs):
......@@ -54,7 +54,7 @@ def main():
plot_comp=NoPlot(),
driver=EasyScipyOptimizeDriver()
)
# view_model(tf)
# n2(tf)
#print(tf.evaluate({'x': desired[:, 0], 'y': desired[:, 1]}))
print(tf.evaluate({'x': optimal[:, 0], 'y': optimal[:, 1]}, disp=False))
#print(tf.evaluate({'x': initial[:, 0], 'y': initial[:, 1]}))
......
......@@ -107,18 +107,18 @@ def main(obj=True, max_con_on=True):
input_keys=[topfarm.x_key, topfarm.y_key, topfarm.type_key, ('obj', obj)],
n_wt=n_wt,
cost_function=aep_func,
output_key="aep",
output_keys="aep",
output_unit="GWh",
objective=obj,
output_val=np.zeros(n_wt),
income_model=True)
output_vals=np.zeros(n_wt),
maximize=True)
comps = [aep_comp] # AEP component is always in the group
if not obj: # if objective is IRR initiate/add irr_comp
irr_comp = CostModelComponent(
input_keys=[topfarm.type_key, 'aep'],
n_wt=n_wt,
cost_function=irr_func,
output_key="irr",
output_keys="irr",
output_unit="%",
objective=True)
comps.append(irr_comp)
......@@ -137,7 +137,7 @@ def main(obj=True, max_con_on=True):
ext_vars=ext_vars)
cost, state, rec = tf.optimize()
# view_model(problem, outfile='ex5_n2.html', show_browser=False)
# n2(problem, outfile='ex5_n2.html', show_browser=False)
end = time.time()
print(end - start)
# %%
......
......@@ -71,18 +71,18 @@ def main():
input_keys=['x', 'y'],
n_wt=n_wt,
cost_function=aep_func,
output_key="aep",
output_keys="aep",
output_unit="GWh",
objective=False,
output_val=np.zeros(n_wt))
output_vals=np.zeros(n_wt))
irr_comp = CostModelComponent(
input_keys=['aep'],
n_wt=n_wt,
cost_function=irr_func,
output_key="irr",
output_keys="irr",
output_unit="%",
objective=True,
income_model=True)
maximize=True)
group = TopFarmGroup([aep_comp, irr_comp])
problem = TopFarmProblem(
design_vars=dict(zip('xy', site.initial_position.T)),
......
......@@ -29,7 +29,7 @@ setup(name='topfarm',
'matplotlib', # for plotting
'numpy', # for numerical calculations
'numpy-financial', # for irr calculations
'openmdao==2.6', # for optimization
'openmdao==3.16', # for optimization
'pytest', # for testing
'pycodestyle', # for testing
'pytest-cov', # for calculating coverage
......
......@@ -26,6 +26,7 @@ from openmdao.api import Problem, IndepVarComp, Group, ParallelGroup,\
ExplicitComponent, ListGenerator, DOEDriver, SimpleGADriver
from openmdao.drivers.doe_generators import DOEGenerator
from openmdao.utils import mpi
from openmdao.core.constants import _SetupStatus
import topfarm
from topfarm.recorders import NestedTopFarmListRecorder,\
TopFarmListRecorder, split_record_id
......@@ -235,21 +236,6 @@ class TopFarmProblem(Problem):
if cost_comp:
self.model.add_subsystem('cost_comp', cost_comp, promotes=['*'])
if expected_cost is None:
expected_cost = self.evaluate()[0]
self._setup_status = 0
if isinstance(driver, EasyDriverBase) and driver.supports_expected_cost is False:
expected_cost = 1
if isinstance(cost_comp, Group) and approx_totals:
cost_comp.approx_totals(**approx_totals)
# Use the assembled Jacobian.
if 'assembled_jac_type' in self.model.cost_comp.options:
self.model.cost_comp.options['assembled_jac_type'] = 'dense'
self.model.cost_comp.linear_solver.assemble_jac = True
else:
self.indeps.add_output('cost')
if len(post_constraints) > 0:
if constraints_as_penalty:
penalty_comp = PostPenaltyComponent(post_constraints, constraints_as_penalty)
......@@ -279,7 +265,23 @@ class TopFarmProblem(Problem):
aggr_comp = AggregatedCost(constraints_as_penalty, constraints, post_constraints)
self.model.add_subsystem('aggr_comp', aggr_comp, promotes=['*'])
# print(expected_cost)
if cost_comp:
if expected_cost is None:
expected_cost = self.evaluate()[0]
if self._metadata:
self._metadata['setup_status'] = 0
if isinstance(driver, EasyDriverBase) and driver.supports_expected_cost is False:
expected_cost = 1
if isinstance(cost_comp, Group) and approx_totals:
cost_comp.approx_totals(**approx_totals)
# Use the assembled Jacobian.
if 'assembled_jac_type' in self.model.cost_comp.options:
self.model.cost_comp.options['assembled_jac_type'] = 'dense'
self.model.cost_comp.linear_solver.assemble_jac = True
else:
self.indeps.add_output('cost')
self.model.add_objective('aggr_cost', scaler=1 / abs(expected_cost))
if plot_comp and not isinstance(plot_comp, NoPlot):
......@@ -341,9 +343,11 @@ class TopFarmProblem(Problem):
topfarm.x_key: x, topfarm.y_key: y, 'c': c}
def setup(self):
if self._setup_status == 0:
if not self._metadata:
Problem.setup(self, check=True)
if self._metadata['setup_status'] == _SetupStatus.PRE_SETUP:
Problem.setup(self, check=True)
if self._setup_status < 2:
if self._metadata['setup_status'] < _SetupStatus.POST_FINAL_SETUP:
with warnings.catch_warnings():
warnings.filterwarnings('error')
try:
......@@ -404,7 +408,7 @@ class TopFarmProblem(Problem):
print("Gradients evaluated in\t%.3fs" % (time.time() - t))
return res
def optimize(self, state={}, disp=False):
def optimize(self, state={}, disp=False, recorder_as_list=False):
"""Run the optimization problem
Parameters
......@@ -415,12 +419,15 @@ class TopFarmProblem(Problem):
The current state is used to unspecified variables
disp : bool, optional
if True, the time used for the optimization is printed
recorder_as_list : bool, optional
if True, returns multiprocessing friendly recorder as list of class and attributes
that can be pickled. Use TopFarmListRecorder().list2recorder to restore TopFarmListRecorder object
Returns
-------
Optimized cost : float
state : dict
recorder : TopFarmListRecorder or NestedTopFarmListRecorder
recorder : TopFarmListRecorder or NestedTopFarmListRecorder or [TopFarmListRecorder.__class__, attributes]
"""
self.load_recorder()
self.update_state(state)
......@@ -454,7 +461,10 @@ class TopFarmProblem(Problem):
best_case_index = int(np.argmin(costs))
best_state = {k: self.recorder[k][best_case_index] for k in self.design_vars}
self.evaluate(best_state)
return self.cost, copy.deepcopy(self.state), self.recorder
if recorder_as_list:
return self.cost, copy.deepcopy(self.state), self.recorder.recorder2list()
else:
return self.cost, copy.deepcopy(self.state), self.recorder
def check_gradients(self, check_all=False, tol=1e-3):
"""Check gradient computations"""
......@@ -535,16 +545,24 @@ class ProblemComponent(ExplicitComponent):
def setup(self):
missing_in_problem_exceptions = ['penalty']
missing_in_problem = (set([c[0] for c in self.parent.indeps._indep_external]) -
set([c[0] for c in self.problem.indeps._indep_external]))
parent_temp = {}
parent_temp.update(self.parent.indeps._static_var_rel2meta)
parent_temp.update(self.parent.indeps._var_rel2meta)
problem_temp = {}
problem_temp.update(self.problem.indeps._static_var_rel2meta)
problem_temp.update(self.problem.indeps._var_rel2meta)
missing_in_problem = (set(parent_temp) -
set(problem_temp))
indepsargs = ['val', 'units']
self.missing_attrs = []
for name, val, kwargs in self.parent.indeps._indep_external:
self.add_input(name, val=val, **{k: kwargs[k] for k in ['units']})
for name, kwargs in parent_temp.items():
kwargs = {k: v for k, v in kwargs.items() if k in indepsargs}
self.add_input(name, **kwargs)
self.missing_attrs.append(name)
if name in missing_in_problem:
if name not in missing_in_problem_exceptions:
self.problem.indeps.add_output(name, val, **kwargs)
self.problem._setup_status = 0 # redo initial setup
self.problem.indeps.add_output(name, **kwargs)
self.problem._setup_status = 1 # 1 -- The `setup` method has been called, but vectors not initialized.
self.problem.setup()
self.add_output('cost', val=0.0)
......@@ -578,7 +596,7 @@ def main():
from topfarm.cost_models.dummy import DummyCost, DummyCostPlotComp
from topfarm.constraint_components.spacing import SpacingConstraint
from topfarm.constraint_components.boundary import XYBoundaryConstraint
from openmdao.api import view_model
from openmdao.api import n2
initial = np.array([[6, 0], [6, -8], [1, 1]]) # initial turbine layouts
optimal = np.array([[2.5, -3], [6, -7], [4.5, -3]]) # optimal turbine layouts
......
......@@ -173,8 +173,14 @@ class BoundaryBaseComp(ConstraintComponent):
def plot(self, ax):
"""Plot boundary"""
ax.plot(self.xy_boundary[:, 0].tolist() + [self.xy_boundary[0, 0]],
self.xy_boundary[:, 1].tolist() + [self.xy_boundary[0, 1]], 'k')
if isinstance(self, MultiPolygonBoundaryComp):
colors = ['--k', 'k']
for bound, io in self.xy_multi_boundary:
ax.plot(np.asarray(bound)[:, 0].tolist() + [np.asarray(bound)[0, 0]],
np.asarray(bound)[:, 1].tolist() + [np.asarray(bound)[0, 1]], colors[io])
else:
ax.plot(self.xy_boundary[:, 0].tolist() + [self.xy_boundary[0, 0]],
self.xy_boundary[:, 1].tolist() + [self.xy_boundary[0, 1]], 'k')
class ConvexBoundaryComp(BoundaryBaseComp):
......@@ -300,7 +306,7 @@ class ConvexBoundaryComp(BoundaryBaseComp):
return self.dfaceDistance_dx, self.dfaceDistance_dy
def satisfy(self, state, pad=1.1):
x, y = [np.asarray(state[xyz], dtype=np.float) for xyz in [topfarm.x_key, topfarm.y_key]]
x, y = [np.asarray(state[xyz], dtype=float) for xyz in [topfarm.x_key, topfarm.y_key]]
dist = self.distances(x, y)
dist = np.where(dist < 0, np.minimum(dist, -.01), dist)
dx, dy = self.gradients(x, y) # independent of position
......@@ -450,7 +456,7 @@ class PolygonBoundaryComp(BoundaryBaseComp):
return np.diagflat(dx), np.diagflat(dy)
def satisfy(self, state, pad=1.1):
x, y = [np.asarray(state[xy], dtype=np.float) for xy in [topfarm.x_key, topfarm.y_key]]
x, y = [np.asarray(state[xy], dtype=float) for xy in [topfarm.x_key, topfarm.y_key]]
dist = self.distances(x, y)
dx, dy = map(np.diag, self.gradients(x, y))
m = dist < 0
......@@ -714,54 +720,54 @@ class MultiPolygonBoundaryComp(PolygonBoundaryComp):
def main():
import matplotlib.pyplot as plt
plt.close('all')
i1 = np.array([[2, 17], [6, 23], [16, 23], [26, 15], [19, 0], [14, 4], [4, 4]])
e1 = np.array([[0, 10], [20, 21], [22, 12], [10, 12], [9, 6], [2, 7]])
i2 = np.array([[12, 13], [14, 17], [18, 15], [17, 10], [15, 11]])
e2 = np.array([[5, 17], [5, 18], [8, 19], [8, 18]])
i3 = np.array([[5, 0], [5, 1], [10, 3], [10, 0]])
e3 = np.array([[6, -1], [6, 18], [7, 18], [7, -1]])
e4 = np.array([[15, 9], [15, 11], [20, 11], [20, 9]])
multi_boundary = [(i1, 'i'), (e1, 'e'), (i2, 'i'), (e2, 'e'), (i3, 'i'), (e3, 'e'), (e4, 'e')]
N_points = 50
xs = np.linspace(-1, 30, N_points)
ys = np.linspace(-1, 30, N_points)
y_grid, x_grid = np.meshgrid(xs, ys)
x = x_grid.ravel()
y = y_grid.ravel()
n_wt = len(x)
MPBC = MultiPolygonBoundaryComp(n_wt, multi_boundary)
distances = MPBC.distances(x, y)
delta = 1e-9
distances2 = MPBC.distances(x + delta, y)
dx_fd = (distances2 - distances) / delta
dx = np.diag(MPBC.gradients(x + delta / 2, y)[0])
plt.figure()
plt.plot(dx_fd, dx, '.')
plt.figure()
for n, bound in enumerate(MPBC.boundaries):
x_bound, y_bound = bound[0].T
x_bound = np.append(x_bound, x_bound[0])
y_bound = np.append(y_bound, y_bound[0])
line, = plt.plot(x_bound, y_bound, label=f'{n}')
plt.plot(x_bound[0], y_bound[0], color=line.get_color(), marker='o')
plt.legend()
plt.grid()
plt.axis('square')
plt.contourf(x_grid, y_grid, distances.reshape(N_points, N_points), 50)
plt.colorbar()
plt.figure()
ax = plt.axes(projection='3d')
ax.contour3D(x.reshape(N_points, N_points), y.reshape(N_points, N_points), distances.reshape(N_points, N_points), 50)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
if __name__ == '__main__':
main()
if __name__ == '__main__':
import matplotlib.pyplot as plt
plt.close('all')
i1 = np.array([[2, 17], [6, 23], [16, 23], [26, 15], [19, 0], [14, 4], [4, 4]])
e1 = np.array([[0, 10], [20, 21], [22, 12], [10, 12], [9, 6], [2, 7]])
i2 = np.array([[12, 13], [14, 17], [18, 15], [17, 10], [15, 11]])
e2 = np.array([[5, 17], [5, 18], [8, 19], [8, 18]])
i3 = np.array([[5, 0], [5, 1], [10, 3], [10, 0]])
e3 = np.array([[6, -1], [6, 18], [7, 18], [7, -1]])
e4 = np.array([[15, 9], [15, 11], [20, 11], [20, 9]])
multi_boundary = [(i1, 'i'), (e1, 'e'), (i2, 'i'), (e2, 'e'), (i3, 'i'), (e3, 'e'), (e4, 'e')]
N_points = 50
xs = np.linspace(-1, 30, N_points)
ys = np.linspace(-1, 30, N_points)
y_grid, x_grid = np.meshgrid(xs, ys)
x = x_grid.ravel()
y = y_grid.ravel()
n_wt = len(x)
MPBC = MultiPolygonBoundaryComp(n_wt, multi_boundary)
distances = MPBC.distances(x, y)
delta = 1e-9
distances2 = MPBC.distances(x + delta, y)
dx_fd = (distances2 - distances) / delta
dx = np.diag(MPBC.gradients(x + delta / 2, y)[0])
plt.figure()
plt.plot(dx_fd, dx, '.')
plt.figure()
for n, bound in enumerate(MPBC.boundaries):
x_bound, y_bound = bound[0].T
x_bound = np.append(x_bound, x_bound[0])
y_bound = np.append(y_bound, y_bound[0])
line, = plt.plot(x_bound, y_bound, label=f'{n}')
plt.plot(x_bound[0], y_bound[0], color=line.get_color(), marker='o')
plt.legend()
plt.grid()
plt.axis('square')
plt.contourf(x_grid, y_grid, distances.reshape(N_points, N_points), 50)
plt.colorbar()
plt.figure()
ax = plt.axes(projection='3d')
ax.contour3D(x.reshape(N_points, N_points), y.reshape(N_points, N_points), distances.reshape(N_points, N_points), 50)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
main()
......@@ -234,7 +234,7 @@ class ConvexBoundaryComp(BoundaryBaseComp):
# partials['boundaryDistances', topfarm.y_key] = dy
def move_inside(self, turbineX, turbineY, turbineZ, pad=1.1):
x, y, z = [np.asarray(xyz, dtype=np.float) for xyz in [turbineX, turbineY, turbineZ]]
x, y, z = [np.asarray(xyz, dtype=float) for xyz in [turbineX, turbineY, turbineZ]]
dist = self.distances(turbineX, turbineY)
dx, dy = self.gradients(x, y) # independent of position
dx = dx[:self.nVertices, 0]
......@@ -390,7 +390,7 @@ class PolygonBoundaryComp(BoundaryBaseComp):
return np.diagflat(dx), np.diagflat(dy)
def move_inside(self, turbineX, turbineY, turbineZ, pad=1.1):
x, y, z = [np.asarray(xyz, dtype=np.float) for xyz in [turbineX, turbineY, turbineZ]]
x, y, z = [np.asarray(xyz, dtype=float) for xyz in [turbineX, turbineY, turbineZ]]
dist = self.distances(turbineX, turbineY)
dx, dy = map(np.diag, self.gradients(x, y))
m = dist < 0
......
......@@ -49,7 +49,7 @@ class CapacityComp(ConstraintComponent):
self.const_id = const_id
def setup(self):
self.add_input(topfarm.type_key, np.zeros(self.n_wt, dtype=np.int))
self.add_input(topfarm.type_key, np.zeros(self.n_wt, dtype=int))
self.add_output('penalty_' + self.const_id, val=0.0)
self.add_output('totalcapacity', val=0.0,
desc='wind farm installed capacity')
......
......@@ -128,7 +128,7 @@ class SpacingComp(ConstraintComponent):
ax.add_artist(circle)
def satisfy(self, state, n_iter=100, step_size=0.1):
x, y = [state[xy].astype(np.float) for xy in [topfarm.x_key, topfarm.y_key]]
x, y = [state[xy].astype(float) for xy in [topfarm.x_key, topfarm.y_key]]
pair_i, pair_j = np.triu_indices(len(x), 1)
for _ in range(n_iter):
dist = self._compute(x, y)
......
......@@ -43,7 +43,7 @@ class CostModelComponent(ExplicitComponent):
maximize : boolean
If True: objective is maximised during optimization\n
If False: Objective is minimized during optimization
output_val : float or array_like
output_vals : float or array_like
Format of output
input_units : list of str
Units of the respective input_keys
......@@ -58,14 +58,10 @@ class CostModelComponent(ExplicitComponent):
if 'output_key' in kwargs:
warnings.warn("""output_key is deprecated; use keyword output_keys instead""",
DeprecationWarning, stacklevel=2)
self.output_key = kwargs['output_key']
# if self.output_key not in output_keys:
output_keys = self.output_key
output_keys = [kwargs['output_key']]
kwargs.pop('output_key')
else:
self.output_key = output_keys[0]
if isinstance(self.output_key, tuple):
self.output_key = self.output_key[0]