Skip to content
Snippets Groups Projects
Commit 0fa0746b authored by Rlahuerta's avatar Rlahuerta
Browse files

ks function update, some fixing into the recorder for IPOPT wrapper

parent ec6828a1
No related branches found
No related tags found
No related merge requests found
Pipeline #35231 failed
Stage: test
Hide whitespace changes
Inline Side-by-side
examples/optimization_course/exclusion_zones.py
+ 2
− 2
View file @ 0fa0746b
......@@ -36,8 +36,8 @@ maximum_water_depth = -52
n_wt = x_init.size
# maxiter = 100
maxiter = 30
# maxiter = 10
# maxiter = 30
maxiter = 10
def aep_func(x: np.ndarray, y: np.ndarray, **kwargs):
......
topfarm/easy_drivers.py
+ 74
− 53
View file @ 0fa0746b
# -*- coding: utf-8 -*-
"""
TODO: docstring
"""
# import copy
import sys
import warnings
from typing import Tuple
# from packaging.version import Version
......@@ -24,7 +28,13 @@ from openmdao.drivers.scipy_optimizer import ScipyOptimizeDriver
from topfarm.drivers.genetic_algorithm_driver import SimpleGADriver
from topfarm.drivers.random_search_driver import RandomSearchDriver
import warnings
try:
from ipyopt import Problem
except ImportError as import_err:
Problem = None
print(f' Warning: {import_err} ')
warnings.filterwarnings("ignore")
CITATIONS = """
......@@ -48,6 +58,10 @@ CITATIONS = """
# noinspection PyCompatibility
class EasyDriverBase(Driver):
"""
TODO: missing-function-docstring
"""
expected_cost = 1
max_iter = None
......@@ -95,6 +109,9 @@ class EasyDriverBase(Driver):
self._cons[name]['linear'] = True
def get_constraints(self, state: dict = None):
"""
TODO: missing-function-docstring
"""
# Pass in new inputs
if state is None:
......@@ -140,8 +157,7 @@ class EasyDriverBase(Driver):
else:
upper_ij = meta_i['upper']
if upper_ij > 1.e16:
upper_ij = 1.e16
upper_ij = min(upper_ij, 1e+16)
if isinstance(meta_i['lower'], np.ndarray):
lower_ij = meta_i['lower'][j]
......@@ -204,6 +220,9 @@ class EasyDriverBase(Driver):
class EasyScipyOptimizeDriver(ScipyOptimizeDriver, EasyDriverBase):
"""
TODO: missing-function-docstring
"""
def __init__(self, optimizer: str = 'SLSQP', maxiter: int = 200, tol: float = 1.e-8, disp: bool = True, **kwargs):
"""
......@@ -325,13 +344,11 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
Finite difference method for approximating the first-order derivatives of the objective function
"""
try:
from ipyopt import Problem
if Problem is not None:
# Store the IpOptWrapper Object
self.ipopt_nlp = Problem
except ImportError:
else:
raise ImportError("""Cannot import ipopt wrapper. Please install ipyopt, e.g. via conda
Windows: no candidate library for now, Linux/OSX: conda install -c conda-forge ipyopt""")
......@@ -347,14 +364,13 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
self._openmdao = False
# Minimum Spacing Constraint Value
self.min_spacing: float = None
self.min_spacing = None
self._list_map_dist: [np.ndarray] = []
self._dist_max = []
self._dist_bounds = []
# TODO: implementation of adaptative ks penalization
self._dist_rho = 40.
# self._dist_rho = 20.
super().__init__(**kwargs)
......@@ -363,10 +379,10 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
self.cite = CITATIONS
self.nvars = 0 # Number of Design Variables
self.xi = np.empty([], dtype=float) # Design Variables Array
self.nvars = 0 # Number of Design Variables
self.xi = np.empty([], dtype=float) # Design Variables Array
self.n_wt = 0 # Number of Turbines
self.n_wt = 0 # Number of Turbines
# Design Variable Array Index
self.state_index = dict(x=np.empty([], dtype=int), y=np.empty([], dtype=int))
......@@ -392,9 +408,9 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
# Delta differences between Upper and Lower bounds Constraints, used for multistart
self.delta_ineq = np.empty([], dtype=float)
self.cst_num = 0 # Number of constraints
self.cst_fun: [{}] = [] # Constraint (function and gradients)
self.cst_spl_fun: [{}] = [] # Constraint (function and gradients) - B-spline or its derivatives
self.cst_num = 0 # Number of constraints
self.cst_fun: [{}] = [] # Constraint (function and gradients)
self.cst_spl_fun: [{}] = [] # Constraint (function and gradients) - B-spline or its derivatives
############################################################
# Full Hessian
......@@ -433,10 +449,6 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
if self.mstart == 1:
self.ipopt_params.update(dict(max_iter=self.maxiter))
def add_min_distance(self, min_spacing: float):
self.min_spacing = min_spacing
def _distance_map(self):
if self.nvars > 0:
......@@ -445,9 +457,6 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
np_arange_idx = np.arange(0, self.n_wt)
for i in range(self.n_wt):
# if i > 0:
# self._list_map_dist.append(np_arange_idx[:i])
idx_i = np.where(i != np_arange_idx)[0]
self._list_map_dist.append(np_arange_idx[idx_i])
......@@ -488,7 +497,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
np_nrm_dist_idx = self._distance_search(np_xi, idx)
np_gx_dist_idx = 1.005 - np_nrm_dist_idx
np_gx_dist_idx = (1. + 0.8 * np.log(np_nrm_dist_idx.shape[0]) / self._dist_rho) - np_nrm_dist_idx
ks_val = self._ks_aggreation(np_gx_dist_idx, stab=self._dist_max[idx], rho=self._dist_rho)
......@@ -498,7 +507,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
return self._fd_diff(self._min_distance_ks_fun, np_xi, idx=idx, dtype='second-order')
def _get_states(self, np_xi: np.ndarray) -> {'x': np.ndarray, 'y': np.ndarray}:
def _get_states(self, np_xi: np.ndarray) -> dict(x=np.ndarray, y=np.ndarray):
return {key_i: np_xi[val_i] for key_i, val_i in self.state_index.items()}
......@@ -531,7 +540,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
self.xi = np.concatenate(list_vars)
self.nvars = self.xi.shape[0]
self.n_wt = np_vars_i.shape[0]
self.n_wt = list_vars[0].shape[0]
#######################################################################
# Box Constraint / lower and upper bounds on the variables xi (Setup)
......@@ -604,7 +613,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
self.opt_settings['disp'] = self.options['disp']
def _bounds_setup(self) -> np.ndarray:
def _bounds_setup(self) -> dict(xi=np.ndarray, lwr=np.ndarray, upp=np.ndarray):
"""Converts bounds to the form required by the solver."""
desvar_vals = self.get_design_var_values()
......@@ -614,7 +623,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
list_bounds_lwr = []
list_bounds_upp = []
for i, (name_i, desvar_i) in enumerate(self._designvars.items()):
for name_i, desvar_i in enumerate(self._designvars.items()):
size_i = desvar_i['global_size'] if desvar_i['distributed'] else desvar_i['size']
np_xi_i = np.zeros(size_i, dtype=float)
......@@ -673,13 +682,13 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
list_cst_lwr = []
list_cst_upp = []
k = 1 # start at 1 since row 0 is the objective. Constraints start at row 1.
lin_i = 0 # counter for linear constraint jacobian
lincons = [] # list of linear constraints
k = 1 # start at 1 since row 0 is the objective. Constraints start at row 1.
lin_i = 0 # counter for linear constraint jacobian
lincons = [] # list of linear constraints
self._obj_and_nlcons = list(self._objs)
for i, (name_i, meta_i) in enumerate(self._cons.items()):
for name_i, meta_i in self._cons.items():
if meta_i['indices'] is not None:
meta_i['size'] = size_i = meta_i['indices'].indexed_src_size
......@@ -743,8 +752,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
else:
upper_ij = meta_i['upper']
if upper_ij > 1.e16:
upper_ij = 1.e16
upper_ij = min(upper_ij, 1.e+16)
if isinstance(meta_i['lower'], np.ndarray):
lower_ij = meta_i['lower'][j]
......@@ -770,7 +778,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
list_cst_lwr += list_cst_lwr_i
list_cst_upp += list_cst_upp_i
list_cst_val_i = [0. for k in list_cst_lwr_i] # unknown g(x) values
list_cst_val_i = [0. for k in list_cst_lwr_i] # unknown g(x) values
self.cst_gxi_bounds[name_i] = pd.DataFrame(dict(lwr=list_cst_lwr_i, val=list_cst_val_i,
upp=list_cst_upp_i, type=list_cst_type_i))
......@@ -801,10 +809,10 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
self.gj_upp, # upper bounds on constraints
self._cst_sparsity_indices(), # number of nonzeros in the constraint Jacobian
self.hessian_sparsity_indices, # number of nonzeros in the Hessian
self._ipopt_objfunc, # objective function
self._ipopt_gradfunc, # gradient of the objective function
self._ipopt_confunc, # constraint function
self._ipopt_confunc_jac # gradient of the constraint function
self._ipopt_objfunc, # objective function
self._ipopt_gradfunc, # gradient of the objective function
self._ipopt_confunc, # constraint function
self._ipopt_confunc_jac # gradient of the constraint function
)
nlp.set(**self.ipopt_params)
......@@ -824,7 +832,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
else:
nsearch = self.mstart
for i in range(nsearch):
for _ in range(nsearch):
np_rd_i = delta_xi * np.random.uniform(size=self.nvars)
np_xi_i = self.perturbation * np_rd_i + x_init
......@@ -922,7 +930,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
self._problem().comm.Bcast(xi_in, root=0)
# Assign Design Variables Values
for i, (name_i, meta_i) in enumerate(self._designvars.items()):
for i, name_i in enumerate(self._designvars.keys()):
self.set_design_var(name_i, dict_states[i])
self._problem().model.run_solve_nonlinear()
......@@ -1023,7 +1031,12 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
Returns
-------
Array with Gradient of the function
Array with Gradient of the function
References
----------
https://www.dam.brown.edu/people/alcyew/handouts/numdiff.pdf
https://math.stackexchange.com/questions/2297749/approximation-of-the-first-derivative-by-writing-taylor-expansions
"""
# FD perturbation value
......@@ -1032,7 +1045,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
# Design Variables Perturbation
np_dx = dtp * xi_new
if dtype == 'backward' or dtype == 'forward':
if dtype in ('backward', 'forward'):
ndiff = 1
else:
ndiff = 2
......@@ -1056,7 +1069,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
else:
fval = fun(xi_new, idx)
for i, x_i in enumerate(xi_new):
for i in range(xi_new.shape[0]):
if dtype == 'backward':
np_xi_i = xi_new.copy()
......@@ -1111,6 +1124,9 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
df_dx_i = ((-1. / 3.) * np_fval_i[0] + (1. / 4.) * fval + (1. / 12.) * np_fval_i[1]) / np_dx[i]
else:
raise NotImplementedError
np_obj_val[i, :] = np_fval_i
np_dobj_val[i] = df_dx_i
......@@ -1256,9 +1272,6 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
return np_cst_val
else:
return
def _ipopt_confunc_jac(self, x_new: np.ndarray, np_jac_out: np.ndarray = None) -> np.ndarray:
"""
Return the value of the IpOpt gradient functions constraints (jacobian matrix).
......@@ -1292,14 +1305,24 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
# self._dist_rho += 0.25
##############################################################################
# Pass in new inputs
# Recorder - openmdao
i = 0
for name, meta in self._designvars.items():
size = meta['size']
self.set_design_var(name, self.xi[i:i + size])
i += size
self._con_cache = self.get_constraint_values()
with RecordingDebugging(self._get_name(), self.iter_count, self) as rec:
self.iter_count += 1
self.topfarm_problem.model.run_solve_nonlinear()
##############################################################################
list_gx_val = []
for i, cst_i in enumerate(self.cst_fun):
for cst_i in self.cst_fun:
if cst_i.get('args') is not None:
name_i = cst_i['args'][0]
idx_i = cst_i['args'][2]
......@@ -1331,9 +1354,6 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
return np_jac_out
else:
return
def run(self) -> bool:
"""
Optimize the problem using selected Scipy optimizer.
......@@ -1360,7 +1380,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
nlp = self._ipopt_nlp_setup()
list_opt = []
for i, np_xi_i in enumerate(np_xi_trials[:self.mstart]):
for np_xi_i in np_xi_trials[:self.mstart]:
list_opt.append(self._solve(nlp, np_xi_i))
if self.mstart == 1:
......@@ -1413,6 +1433,7 @@ class EasyIpOptOptimizeDriver(drv.ScipyOptimizeDriver, EasyDriverBase):
def _get_name(self):
"""Override to add str"""
return "Optimize_IpOpt"
......@@ -1461,7 +1482,7 @@ except ModuleNotFoundError:
class EasySimpleGADriver(SimpleGADriver, EasyDriverBase):
def __init__(self, max_gen: int = 100, pop_size: int = 25, Pm: float = None, Pc=.5, elitism=True,
def __init__(self, max_gen: int = 100, pop_size: int = 25, Pm: float = None, Pc: float = .5, elitism: bool = True,
bits={}, debug_print=[], run_parallel=False, random_state=None):
"""SimpleGA driver with optional arguments
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment