diff --git a/docs/wisdem/nrelcsm/theory.rst b/docs/wisdem/nrelcsm/theory.rst index 64d5b5a23..eb0d37f85 100644 --- a/docs/wisdem/nrelcsm/theory.rst +++ b/docs/wisdem/nrelcsm/theory.rst @@ -42,12 +42,12 @@ To obtain the blade mass in kilograms and cost in USD from the rotor diameter in b &= (see below)\\ k_c &= 14.6 -Where :math:`D_{rotor}` is the rotor diameter and :math:`k_c` is determined by: +Where :math:`D_{rotor}` is the rotor diameter and :math:`b` is determined by: -* If turbine class I and blade DOES have carbon fiber spar caps, :math:`k_c=2.47` -* If turbine class I and blade DOES NOT have carbon fiber spar caps, :math:`k_c=2.54` -* If turbine class II+ and blade DOES have carbon fiber spar caps, :math:`k_c=2.44` -* If turbine class II+ and blade DOES NOT have carbon fiber spar caps, :math:`k_c=2.50` +* If turbine class I and blade DOES have carbon fiber spar caps, :math:`b=2.47` +* If turbine class I and blade DOES NOT have carbon fiber spar caps, :math:`b=2.54` +* If turbine class II+ and blade DOES have carbon fiber spar caps, :math:`b=2.44` +* If turbine class II+ and blade DOES NOT have carbon fiber spar caps, :math:`b=2.50` * User override of exponent value For variable names access to override the default values see the :ref:`csmsource`. diff --git a/wisdem/commonse/utilities.py b/wisdem/commonse/utilities.py index affb87d1c..1bd1968f5 100644 --- a/wisdem/commonse/utilities.py +++ b/wisdem/commonse/utilities.py @@ -31,7 +31,7 @@ def get_modal_coefficients(x, y, deg=6): return p6 -def get_mode_shapes(r, freqs, xdsp, ydsp, zdsp, xmpf, ympf, zmpf): +def get_xy_mode_shapes(r, freqs, xdsp, ydsp, zdsp, xmpf, ympf, zmpf): # Number of frequencies and modes nfreq = len(freqs) @@ -61,8 +61,6 @@ def get_mode_shapes(r, freqs, xdsp, ydsp, zdsp, xmpf, ympf, zmpf): mshapes_y[iy,:] = ypolys[m,:] freq_y[ iy ] = freqs[m] iy += 1 - else: - print('Warning: Unknown mode shape') return freq_x, freq_y, mshapes_x, mshapes_y diff --git a/wisdem/commonse/vertical_cylinder.py b/wisdem/commonse/vertical_cylinder.py index 703a43a7a..6606ca3f2 100644 --- a/wisdem/commonse/vertical_cylinder.py +++ b/wisdem/commonse/vertical_cylinder.py @@ -545,7 +545,8 @@ def compute(self, inputs, outputs): Mmethod = 1 lump = 0 shift = 0.0 - cylinder.enableDynamics(NFREQ, Mmethod, lump, frame3dd_opt['tol'], shift) + # Run twice the number of modes to ensure that we can ignore the torsional modes and still get the desired number of fore-aft, side-side modes + cylinder.enableDynamics(2*NFREQ, Mmethod, lump, frame3dd_opt['tol'], shift) # ---------------------------- # ------ static load case 1 ------------ @@ -592,15 +593,16 @@ def compute(self, inputs, outputs): # natural frequncies outputs['f1'] = modal.freq[0] outputs['f2'] = modal.freq[1] - outputs['freqs'] = modal.freq + outputs['freqs'] = modal.freq[:NFREQ] # Get all mode shapes in batch - freq_x, freq_y, mshapes_x, mshapes_y = util.get_mode_shapes(z, modal.freq, modal.xdsp, modal.ydsp, modal.zdsp, modal.xmpf, modal.ympf, modal.zmpf) - outputs['x_mode_freqs'] = freq_x - outputs['y_mode_freqs'] = freq_y - outputs['x_mode_shapes'] = mshapes_x - outputs['y_mode_shapes'] = mshapes_y - + NFREQ2 = int(NFREQ/2) + freq_x, freq_y, mshapes_x, mshapes_y = util.get_xy_mode_shapes(z, modal.freq, modal.xdsp, modal.ydsp, modal.zdsp, modal.xmpf, modal.ympf, modal.zmpf) + outputs['x_mode_freqs'] = freq_x[:NFREQ2] + outputs['y_mode_freqs'] = freq_y[:NFREQ2] + outputs['x_mode_shapes'] = mshapes_x[:NFREQ2,:] + outputs['y_mode_shapes'] = mshapes_y[:NFREQ2,:] + # deflections due to loading (from cylinder top and wind/wave loads) outputs['top_deflection'] = displacements.dx[iCase, n-1] # in yaw-aligned direction diff --git a/wisdem/glue_code/gc_PoseOptimization.py b/wisdem/glue_code/gc_PoseOptimization.py new file mode 100644 index 000000000..716556056 --- /dev/null +++ b/wisdem/glue_code/gc_PoseOptimization.py @@ -0,0 +1,392 @@ +import numpy as np +import openmdao.api as om +import os + +class PoseOptimization(object): + def __init__(self, modeling_options, analysis_options): + self.modeling = modeling_options + self.opt = analysis_options + self.blade_opt = self.opt['optimization_variables']['blade'] + self.tower_opt = self.opt['optimization_variables']['tower'] + self.control_opt = self.opt['optimization_variables']['control'] + + + def get_number_design_variables(self): + # Determine the number of design variables + n_DV = 0 + + if self.blade_opt['aero_shape']['twist']['flag']: + n_DV += self.blade_opt['aero_shape']['twist']['n_opt'] - 2 + if self.blade_opt['aero_shape']['chord']['flag']: + n_DV += self.blade_opt['aero_shape']['chord']['n_opt'] - 3 + if self.blade_opt['aero_shape']['af_positions']['flag']: + n_DV += self.modeling['blade']['n_af_span'] - self.blade_opt['aero_shape']['af_positions']['af_start'] - 1 + if self.blade_opt['structure']['spar_cap_ss']['flag']: + n_DV += self.blade_opt['structure']['spar_cap_ss']['n_opt'] - 2 + if self.blade_opt['structure']['spar_cap_ps']['flag'] and not self.blade_opt['structure']['spar_cap_ps']['equal_to_suction']: + n_DV += self.blade_opt['structure']['spar_cap_ps']['n_opt'] - 2 + if self.opt['optimization_variables']['control']['tsr']['flag']: + n_DV += 1 + if self.opt['optimization_variables']['control']['servo']['pitch_control']['flag']: + n_DV += 2 + if self.opt['optimization_variables']['control']['servo']['torque_control']['flag']: + n_DV += 2 + if self.tower_opt['outer_diameter']['flag']: + n_DV += self.modeling['tower']['n_height'] + if self.tower_opt['layer_thickness']['flag']: + n_DV += (self.modeling['tower']['n_height'] - 1) * self.modeling['tower']['n_layers'] + + if self.opt['driver']['form'] == 'central': + n_DV *= 2 + + return n_DV + + + def _get_step_size(self): + # If a step size for the driver-level finite differencing is provided, use that step size. Otherwise use a default value. + return (1.e-6 if not 'step_size' in self.opt['driver'] else self.opt['driver']['step_size']) + + + def set_driver(self, wt_opt): + folder_output = self.opt['general']['folder_output'] + + step_size = self._get_step_size() + + # Solver has specific meaning in OpenMDAO + wt_opt.model.approx_totals(method='fd', step=step_size, form=self.opt['driver']['form']) + + # Set optimization solver and options. First, Scipy's SLSQP + if self.opt['driver']['solver'] == 'SLSQP': + wt_opt.driver = om.ScipyOptimizeDriver() + wt_opt.driver.options['optimizer'] = self.opt['driver']['solver'] + wt_opt.driver.options['tol'] = self.opt['driver']['tol'] + wt_opt.driver.options['maxiter'] = self.opt['driver']['max_iter'] + + # The next two optimization methods require pyOptSparse. + elif self.opt['driver']['solver'] == 'CONMIN': + try: + from openmdao.api import pyOptSparseDriver + except: + exit('You requested the optimization solver CONMIN, but you have not installed the pyOptSparseDriver. Please do so and rerun.') + wt_opt.driver = pyOptSparseDriver() + wt_opt.driver.options['optimizer'] = self.opt['driver']['solver'] + wt_opt.driver.opt_settings['ITMAX']= self.opt['driver']['max_iter'] + + elif self.opt['driver']['solver'] == 'SNOPT': + try: + from openmdao.api import pyOptSparseDriver + except: + exit('You requested the optimization solver SNOPT, but you have not installed the pyOptSparseDriver. Please do so and rerun.') + wt_opt.driver = pyOptSparseDriver() + try: + wt_opt.driver.options['optimizer'] = self.opt['driver']['solver'] + except: + exit('You requested the optimization solver SNOPT, but you have not installed it within the pyOptSparseDriver. Please do so and rerun.') + wt_opt.driver.opt_settings['Major optimality tolerance'] = float(self.opt['driver']['tol']) + wt_opt.driver.opt_settings['Major iterations limit'] = int(self.opt['driver']['max_major_iter']) + wt_opt.driver.opt_settings['Iterations limit'] = int(self.opt['driver']['max_minor_iter']) + wt_opt.driver.opt_settings['Major feasibility tolerance'] = float(self.opt['driver']['tol']) + wt_opt.driver.opt_settings['Summary file'] = os.path.join(folder_output, 'SNOPT_Summary_file.txt') + wt_opt.driver.opt_settings['Print file'] = os.path.join(folder_output, 'SNOPT_Print_file.txt') + if 'hist_file_name' in self.opt['driver']: + wt_opt.driver.hist_file = self.opt['driver']['hist_file_name'] + if 'verify_level' in self.opt['driver']: + wt_opt.driver.opt_settings['Verify level'] = self.opt['driver']['verify_level'] + # wt_opt.driver.declare_coloring() + if 'hotstart_file' in self.opt['driver']: + wt_opt.driver.hotstart_file = self.opt['driver']['hotstart_file'] + + else: + exit('The optimizer ' + self.opt['driver']['solver'] + 'is not yet supported!') + + return wt_opt + + + def set_objective(self, wt_opt): + + # Set merit figure. Each objective has its own scaling. + if self.opt['merit_figure'] == 'AEP': + wt_opt.model.add_objective('sse.AEP', ref = -1.e6) + elif self.opt['merit_figure'] == 'blade_mass': + wt_opt.model.add_objective('elastic.precomp.blade_mass', ref = 1.e4) + elif self.opt['merit_figure'] == 'LCOE': + wt_opt.model.add_objective('financese.lcoe', ref = 0.1) + elif self.opt['merit_figure'] == 'blade_tip_deflection': + wt_opt.model.add_objective('tcons.tip_deflection_ratio') + elif self.opt['merit_figure'] == 'tower_mass': + wt_opt.model.add_objective('towerse.tower_mass') + elif self.opt['merit_figure'] == 'tower_cost': + wt_opt.model.add_objective('tcc.tower_cost') + elif self.opt['merit_figure'] == 'Cp': + if self.modeling['Analysis_Flags']['ServoSE']: + wt_opt.model.add_objective('sse.powercurve.Cp_regII', ref = -1.) + else: + wt_opt.model.add_objective('ccblade.CP', ref = -1.) + else: + exit('The merit figure ' + self.opt['merit_figure'] + ' is not supported.') + + return wt_opt + + + def set_design_variables(self, wt_opt, wt_init): + + # Set optimization design variables. + + if self.blade_opt['aero_shape']['twist']['flag']: + indices = range(2, self.blade_opt['aero_shape']['twist']['n_opt']) + wt_opt.model.add_design_var('blade.opt_var.twist_opt_gain', indices = indices, lower=0., upper=1.) + + chord_options = self.blade_opt['aero_shape']['chord'] + if chord_options['flag']: + indices = range(3, chord_options['n_opt'] - 1) + wt_opt.model.add_design_var('blade.opt_var.chord_opt_gain', indices = indices, lower=chord_options['min_gain'], upper=chord_options['max_gain']) + + if self.blade_opt['aero_shape']['af_positions']['flag']: + n_af = self.modeling['blade']['n_af_span'] + indices = range(self.blade_opt['aero_shape']['af_positions']['af_start'],n_af - 1) + af_pos_init = wt_init['components']['blade']['outer_shape_bem']['airfoil_position']['grid'] + step_size = self._get_step_size() + lb_af = np.zeros(n_af) + ub_af = np.zeros(n_af) + for i in range(1,indices[0]): + lb_af[i] = ub_af[i] = af_pos_init[i] + for i in indices: + lb_af[i] = 0.5*(af_pos_init[i-1] + af_pos_init[i]) + step_size + ub_af[i] = 0.5*(af_pos_init[i+1] + af_pos_init[i]) - step_size + lb_af[-1] = ub_af[-1] = 1. + wt_opt.model.add_design_var('blade.opt_var.af_position', indices = indices, lower=lb_af[indices], upper=ub_af[indices]) + + spar_cap_ss_options = self.blade_opt['structure']['spar_cap_ss'] + if spar_cap_ss_options['flag']: + indices = range(1,spar_cap_ss_options['n_opt'] - 1) + wt_opt.model.add_design_var('blade.opt_var.spar_cap_ss_opt_gain', indices = indices, lower=spar_cap_ss_options['min_gain'], upper=spar_cap_ss_options['max_gain']) + + # Only add the pressure side design variables if we do set + # `equal_to_suction` as False in the optimization yaml. + spar_cap_ps_options = self.blade_opt['structure']['spar_cap_ps'] + if spar_cap_ps_options['flag'] and not spar_cap_ps_options['equal_to_suction']: + indices = range(1, spar_cap_ps_options['n_opt'] - 1) + wt_opt.model.add_design_var('blade.opt_var.spar_cap_ps_opt_gain', indices = indices, lower=spar_cap_ps_options['min_gain'], upper=spar_cap_ps_options['max_gain']) + + if self.tower_opt['outer_diameter']['flag']: + wt_opt.model.add_design_var('tower.diameter', lower=self.tower_opt['outer_diameter']['lower_bound'], upper=self.tower_opt['outer_diameter']['upper_bound'], ref=5.) + + if self.tower_opt['layer_thickness']['flag']: + wt_opt.model.add_design_var('tower.layer_thickness', lower=self.tower_opt['layer_thickness']['lower_bound'], upper=self.tower_opt['layer_thickness']['upper_bound'], ref=1e-2) + + # -- Control -- + if self.control_opt['tsr']['flag']: + wt_opt.model.add_design_var('opt_var.tsr_opt_gain', lower=self.control_opt['tsr']['min_gain'], + upper=self.control_opt['tsr']['max_gain']) + if self.control_opt['servo']['pitch_control']['flag']: + wt_opt.model.add_design_var('control.PC_omega', lower=self.control_opt['servo']['pitch_control']['omega_min'], + upper=self.control_opt['servo']['pitch_control']['omega_max']) + wt_opt.model.add_design_var('control.PC_zeta', lower=self.control_opt['servo']['pitch_control']['zeta_min'], + upper=self.control_opt['servo']['pitch_control']['zeta_max']) + if self.control_opt['servo']['torque_control']['flag']: + wt_opt.model.add_design_var('control.VS_omega', lower=self.control_opt['servo']['torque_control']['omega_min'], + upper=self.control_opt['servo']['torque_control']['omega_max']) + wt_opt.model.add_design_var('control.VS_zeta', lower=self.control_opt['servo']['torque_control']['zeta_min'], + upper=self.control_opt['servo']['torque_control']['zeta_max']) + + return wt_opt + + + def set_constraints(self, wt_opt): + + # Set non-linear constraints + blade_constraints = self.opt['constraints']['blade'] + if blade_constraints['strains_spar_cap_ss']['flag']: + if self.blade_opt['structure']['spar_cap_ss']['flag']: + wt_opt.model.add_constraint('rlds.constr.constr_max_strainU_spar', upper= 1.0) + else: + print('WARNING: the strains of the suction-side spar cap are set to be constrained, but spar cap thickness is not an active design variable. The constraint is not enforced.') + + if blade_constraints['strains_spar_cap_ps']['flag']: + if self.blade_opt['structure']['spar_cap_ps']['flag'] or self.blade_opt['structure']['spar_cap_ps']['equal_to_suction']: + wt_opt.model.add_constraint('rlds.constr.constr_max_strainL_spar', upper= 1.0) + else: + print('WARNING: the strains of the pressure-side spar cap are set to be constrained, but spar cap thickness is not an active design variable. The constraint is not enforced.') + + if blade_constraints['stall']['flag']: + if self.blade_opt['aero_shape']['twist']['flag']: + wt_opt.model.add_constraint('stall_check.no_stall_constraint', upper= 1.0) + else: + print('WARNING: the margin to stall is set to be constrained, but twist is not an active design variable. The constraint is not enforced.') + + if blade_constraints['tip_deflection']['flag']: + if self.blade_opt['structure']['spar_cap_ss']['flag'] or self.blade_opt['structure']['spar_cap_ps']['flag']: + wt_opt.model.add_constraint('tcons.tip_deflection_ratio', upper= blade_constraints['tip_deflection']['ratio']) + else: + print('WARNING: the tip deflection is set to be constrained, but spar caps thickness is not an active design variable. The constraint is not enforced.') + + if blade_constraints['chord']['flag']: + if self.blade_opt['aero_shape']['chord']['flag']: + wt_opt.model.add_constraint('blade.pa.max_chord_constr', upper= 1.0) + else: + print('WARNING: the max chord is set to be constrained, but chord is not an active design variable. The constraint is not enforced.') + + if blade_constraints['frequency']['flap_above_3P']: + if self.blade_opt['structure']['spar_cap_ss']['flag'] or self.blade_opt['structure']['spar_cap_ps']['flag']: + wt_opt.model.add_constraint('rlds.constr.constr_flap_f_margin', upper= 0.0) + else: + print('WARNING: the blade flap frequencies are set to be constrained, but spar caps thickness is not an active design variable. The constraint is not enforced.') + + if blade_constraints['frequency']['edge_above_3P']: + wt_opt.model.add_constraint('rlds.constr.constr_edge_f_margin', upper= 0.0) + + if blade_constraints['rail_transport']['flag']: + if blade_constraints['rail_transport']['8_axle']: + wt_opt.model.add_constraint('elastic.rail.constr_LV_8axle_horiz', lower = 0.8, upper= 1.0) + wt_opt.model.add_constraint('elastic.rail.constr_strainPS', upper= 1.0) + wt_opt.model.add_constraint('elastic.rail.constr_strainSS', upper= 1.0) + elif blade_constraints['rail_transport']['4_axle']: + wt_opt.model.add_constraint('elastic.rail.constr_LV_4axle_horiz', upper= 1.0) + else: + exit('You have activated the rail transport constraint module. Please define whether you want to model 4- or 8-axle flatcars.') + + if self.opt['constraints']['blade']['moment_coefficient']['flag']: + wt_opt.model.add_constraint('ccblade.CM', lower= self.opt['constraints']['blade']['moment_coefficient']['min'], upper= self.opt['constraints']['blade']['moment_coefficient']['max']) + if self.opt['constraints']['blade']['match_cl_cd']['flag_cl'] or self.opt['constraints']['blade']['match_cl_cd']['flag_cd']: + data_target = np.loadtxt(self.opt['constraints']['blade']['match_cl_cd']['filename']) + eta_opt = np.linspace(0., 1., self.opt['optimization_variables']['blade']['aero_shape']['twist']['n_opt']) + target_cl = np.interp(eta_opt, data_target[:,0], data_target[:,3]) + target_cd = np.interp(eta_opt, data_target[:,0], data_target[:,4]) + eps_cl = 1.e-2 + if self.opt['constraints']['blade']['match_cl_cd']['flag_cl']: + wt_opt.model.add_constraint('ccblade.cl_n_opt', lower = target_cl-eps_cl, upper = target_cl+eps_cl) + if self.opt['constraints']['blade']['match_cl_cd']['flag_cd']: + wt_opt.model.add_constraint('ccblade.cd_n_opt', lower = target_cd-eps_cl, upper = target_cd+eps_cl) + if self.opt['constraints']['blade']['match_L_D']['flag_L'] or self.opt['constraints']['blade']['match_L_D']['flag_D']: + data_target = np.loadtxt(self.opt['constraints']['blade']['match_L_D']['filename']) + eta_opt = np.linspace(0., 1., self.opt['optimization_variables']['blade']['aero_shape']['twist']['n_opt']) + target_L = np.interp(eta_opt, data_target[:,0], data_target[:,7]) + target_D = np.interp(eta_opt, data_target[:,0], data_target[:,8]) + eps_L = 1.e+2 + if self.opt['constraints']['blade']['match_L_D']['flag_L']: + wt_opt.model.add_constraint('ccblade.L_n_opt', lower = target_L-eps_L, upper = target_L+eps_L) + if self.opt['constraints']['blade']['match_L_D']['flag_D']: + wt_opt.model.add_constraint('ccblade.D_n_opt', lower = target_D-eps_L, upper = target_D+eps_L) + + tower_constraints = self.opt['constraints']['tower'] + if tower_constraints['height_constraint']['flag']: + wt_opt.model.add_constraint('towerse.height_constraint', + lower=tower_constraints['height_constraint']['lower_bound'], + upper=tower_constraints['height_constraint']['upper_bound']) + + if tower_constraints['stress']['flag']: + wt_opt.model.add_constraint('towerse.post.stress', upper=1.0) + + if tower_constraints['global_buckling']['flag']: + wt_opt.model.add_constraint('towerse.post.global_buckling', upper=1.0) + + if tower_constraints['shell_buckling']['flag']: + wt_opt.model.add_constraint('towerse.post.shell_buckling', upper=1.0) + + if tower_constraints['constr_d_to_t']['flag']: + wt_opt.model.add_constraint('towerse.constr_d_to_t', upper=0.0) + + if tower_constraints['constr_taper']['flag']: + wt_opt.model.add_constraint('towerse.constr_taper', lower=0.0) + + if tower_constraints['slope']['flag']: + wt_opt.model.add_constraint('towerse.slope', upper=1.0) + + if tower_constraints['frequency_1']['flag']: + wt_opt.model.add_constraint('towerse.tower.f1', + lower=tower_constraints['frequency_1']['lower_bound'], + upper=tower_constraints['frequency_1']['upper_bound']) + + #control_constraints = self.opt['constraints']['control'] + return wt_opt + + + def set_recorders(self, wt_opt): + folder_output = self.opt['general']['folder_output'] + + # Set recorder on the OpenMDAO driver level using the `optimization_log` + # filename supplied in the optimization yaml + if self.opt['recorder']['flag']: + recorder = om.SqliteRecorder(os.path.join(folder_output, self.opt['recorder']['file_name'])) + wt_opt.driver.add_recorder(recorder) + wt_opt.add_recorder(recorder) + + wt_opt.driver.recording_options['excludes'] = ['*_df'] + wt_opt.driver.recording_options['record_constraints'] = True + wt_opt.driver.recording_options['record_desvars'] = True + wt_opt.driver.recording_options['record_objectives'] = True + + return wt_opt + + + def set_initial(self, wt_opt, wt_init): + + wt_opt['blade.opt_var.s_opt_twist'] = np.linspace(0., 1., self.blade_opt['aero_shape']['twist']['n_opt']) + if self.blade_opt['aero_shape']['twist']['flag']: + init_twist_opt = np.interp(wt_opt['blade.opt_var.s_opt_twist'], wt_init['components']['blade']['outer_shape_bem']['twist']['grid'], wt_init['components']['blade']['outer_shape_bem']['twist']['values']) + lb_twist = np.array(self.blade_opt['aero_shape']['twist']['lower_bound']) + ub_twist = np.array(self.blade_opt['aero_shape']['twist']['upper_bound']) + wt_opt['blade.opt_var.twist_opt_gain'] = (init_twist_opt - lb_twist) / (ub_twist - lb_twist) + if max(wt_opt['blade.opt_var.twist_opt_gain']) > 1. or min(wt_opt['blade.opt_var.twist_opt_gain']) < 0.: + print('Warning: the initial twist violates the upper or lower bounds of the twist design variables.') + + blade_constraints = self.opt['constraints']['blade'] + wt_opt['blade.opt_var.s_opt_chord'] = np.linspace(0., 1., self.blade_opt['aero_shape']['chord']['n_opt']) + wt_opt['blade.ps.s_opt_spar_cap_ss'] = np.linspace(0., 1., self.blade_opt['structure']['spar_cap_ss']['n_opt']) + wt_opt['blade.ps.s_opt_spar_cap_ps'] = np.linspace(0., 1., self.blade_opt['structure']['spar_cap_ps']['n_opt']) + wt_opt['rlds.constr.max_strainU_spar'] = blade_constraints['strains_spar_cap_ss']['max'] + wt_opt['rlds.constr.max_strainL_spar'] = blade_constraints['strains_spar_cap_ps']['max'] + wt_opt['stall_check.stall_margin'] = blade_constraints['stall']['margin'] * 180. / np.pi + + return wt_opt + + + def set_restart(self, wt_opt): + if 'warmstart_file' in self.opt['driver']: + + # Directly read the pyoptsparse sqlite db file + from pyoptsparse import SqliteDict + db = SqliteDict(self.opt['driver']['warmstart_file']) + + # Grab the last iteration's design variables + last_key = db['last'] + desvars = db[last_key]['xuser'] + + # Obtain the already-setup OM problem's design variables + if wt_opt.model._static_mode: + design_vars = wt_opt.model._static_design_vars + else: + design_vars = wt_opt.model._design_vars + + # Get the absolute names from the promoted names within the OM model. + # We need this because the pyoptsparse db has the absolute names for + # variables but the OM model uses the promoted names. + prom2abs = wt_opt.model._var_allprocs_prom2abs_list['output'] + abs2prom = {} + for key in design_vars: + abs2prom[prom2abs[key][0]] = key + + # Loop through each design variable + for key in desvars: + prom_key = abs2prom[key] + + # Scale each DV based on the OM scaling from the problem. + # This assumes we're running the same problem with the same scaling + scaler = design_vars[prom_key]['scaler'] + adder = design_vars[prom_key]['adder'] + + if scaler is None: + scaler = 1.0 + if adder is None: + adder = 0.0 + + scaled_dv = desvars[key] / scaler - adder + + # Special handling for blade twist as we only have the + # last few control points as design variables + if 'twist_opt_gain' in key: + wt_opt[key][2:] = scaled_dv + else: + wt_opt[key][:] = scaled_dv + + return wt_opt diff --git a/wisdem/glue_code/glue_code.py b/wisdem/glue_code/glue_code.py index 4df63cb21..29fed7a86 100644 --- a/wisdem/glue_code/glue_code.py +++ b/wisdem/glue_code/glue_code.py @@ -295,7 +295,6 @@ def setup(self): self.connect('nacelle.lss_diameter', 'drivese.bear1.D_shaft', src_indices=[0]) self.connect('nacelle.lss_diameter', 'drivese.bear2.D_shaft', src_indices=[-1]) self.connect('nacelle.uptower', 'drivese.uptower') - #self.connect('nacelle.gearbox_efficiency', 'drivese.gearbox_efficiency') self.connect('nacelle.brake_mass_user', 'drivese.brake_mass_user') self.connect('nacelle.hvac_mass_coeff', 'drivese.hvac_mass_coeff') self.connect('nacelle.converter_mass_user', 'drivese.converter_mass_user') diff --git a/wisdem/glue_code/runWISDEM.py b/wisdem/glue_code/runWISDEM.py index 48da8782b..6d7fe1f10 100644 --- a/wisdem/glue_code/runWISDEM.py +++ b/wisdem/glue_code/runWISDEM.py @@ -3,6 +3,7 @@ import openmdao.api as om from wisdem.glue_code.gc_LoadInputs import WindTurbineOntologyPython from wisdem.glue_code.gc_WT_InitModel import yaml2openmdao +from wisdem.glue_code.gc_PoseOptimization import PoseOptimization from wisdem.glue_code.glue_code import WindPark from wisdem.commonse.mpi_tools import MPI from wisdem.commonse import fileIO @@ -24,36 +25,12 @@ def run_wisdem(fname_wt_input, fname_modeling_options, fname_opt_options, overri # Initialize openmdao problem. If running with multiple processors in MPI, use parallel finite differencing equal to the number of cores used. # Otherwise, initialize the WindPark system normally. Get the rank number for parallelization. We only print output files using the root processor. - blade_opt_options = opt_options['optimization_variables']['blade'] - tower_opt_options = opt_options['optimization_variables']['tower'] - control_opt_options = opt_options['optimization_variables']['control'] + myopt = PoseOptimization(modeling_options, opt_options) + if MPI: - # Determine the number of design variables - n_DV = 0 - if blade_opt_options['aero_shape']['twist']['flag']: - n_DV += blade_opt_options['aero_shape']['twist']['n_opt'] - 2 - if blade_opt_options['aero_shape']['chord']['flag']: - n_DV += blade_opt_options['aero_shape']['chord']['n_opt'] - 3 - if blade_opt_options['aero_shape']['af_positions']['flag']: - n_DV += modeling_options['blade']['n_af_span'] - blade_opt_options['aero_shape']['af_positions']['af_start'] - 1 - if blade_opt_options['structure']['spar_cap_ss']['flag']: - n_DV += blade_opt_options['structure']['spar_cap_ss']['n_opt'] - 2 - if blade_opt_options['structure']['spar_cap_ps']['flag'] and not blade_opt_options['structure']['spar_cap_ps']['equal_to_suction']: - n_DV += blade_opt_options['structure']['spar_cap_ps']['n_opt'] - 2 - if opt_options['optimization_variables']['control']['tsr']['flag']: - n_DV += 1 - if opt_options['optimization_variables']['control']['servo']['pitch_control']['flag']: - n_DV += 2 - if opt_options['optimization_variables']['control']['servo']['torque_control']['flag']: - n_DV += 2 - if tower_opt_options['outer_diameter']['flag']: - n_DV += modeling_options['tower']['n_height'] - if tower_opt_options['layer_thickness']['flag']: - n_DV += (modeling_options['tower']['n_height'] - 1) * modeling_options['tower']['n_layers'] - - if opt_options['driver']['form'] == 'central': - n_DV *= 2 + n_DV = myopt.get_number_design_variables() + # Extract the number of cores available max_cores = MPI.COMM_WORLD.Get_size() @@ -88,276 +65,18 @@ def run_wisdem(fname_wt_input, fname_modeling_options, fname_opt_options, overri # If at least one of the design variables is active, setup an optimization if opt_options['opt_flag']: - # If a step size for the driver-level finite differencing is provided, use that step size. Otherwise use a default value. - if 'step_size' in opt_options['driver']: - step_size = opt_options['driver']['step_size'] - else: - step_size = 1.e-6 - - # Solver has specific meaning in OpenMDAO - wt_opt.model.approx_totals(method='fd', step=step_size, form=opt_options['driver']['form']) - - # Set optimization solver and options. First, Scipy's SLSQP - if opt_options['driver']['solver'] == 'SLSQP': - wt_opt.driver = om.ScipyOptimizeDriver() - wt_opt.driver.options['optimizer'] = opt_options['driver']['solver'] - wt_opt.driver.options['tol'] = opt_options['driver']['tol'] - wt_opt.driver.options['maxiter'] = opt_options['driver']['max_iter'] - - # The next two optimization methods require pyOptSparse. - elif opt_options['driver']['solver'] == 'CONMIN': - try: - from openmdao.api import pyOptSparseDriver - except: - exit('You requested the optimization solver CONMIN, but you have not installed the pyOptSparseDriver. Please do so and rerun.') - wt_opt.driver = pyOptSparseDriver() - wt_opt.driver.options['optimizer'] = opt_options['driver']['solver'] - wt_opt.driver.opt_settings['ITMAX']= opt_options['driver']['max_iter'] - - elif opt_options['driver']['solver'] == 'SNOPT': - try: - from openmdao.api import pyOptSparseDriver - except: - exit('You requested the optimization solver SNOPT, but you have not installed the pyOptSparseDriver. Please do so and rerun.') - wt_opt.driver = pyOptSparseDriver() - try: - wt_opt.driver.options['optimizer'] = opt_options['driver']['solver'] - except: - exit('You requested the optimization solver SNOPT, but you have not installed it within the pyOptSparseDriver. Please do so and rerun.') - wt_opt.driver.opt_settings['Major optimality tolerance'] = float(opt_options['driver']['tol']) - wt_opt.driver.opt_settings['Major iterations limit'] = int(opt_options['driver']['max_major_iter']) - wt_opt.driver.opt_settings['Iterations limit'] = int(opt_options['driver']['max_minor_iter']) - wt_opt.driver.opt_settings['Major feasibility tolerance'] = float(opt_options['driver']['tol']) - wt_opt.driver.opt_settings['Summary file'] = os.path.join(folder_output, 'SNOPT_Summary_file.txt') - wt_opt.driver.opt_settings['Print file'] = os.path.join(folder_output, 'SNOPT_Print_file.txt') - if 'hist_file_name' in opt_options['driver']: - wt_opt.driver.hist_file = opt_options['driver']['hist_file_name'] - if 'verify_level' in opt_options['driver']: - wt_opt.driver.opt_settings['Verify level'] = opt_options['driver']['verify_level'] - # wt_opt.driver.declare_coloring() - if 'hotstart_file' in opt_options['driver']: - wt_opt.driver.hotstart_file = opt_options['driver']['hotstart_file'] - - else: - exit('The optimizer ' + opt_options['driver']['solver'] + 'is not yet supported!') - - # Set merit figure. Each objective has its own scaling. - if opt_options['merit_figure'] == 'AEP': - wt_opt.model.add_objective('sse.AEP', ref = -1.e6) - elif opt_options['merit_figure'] == 'blade_mass': - wt_opt.model.add_objective('elastic.precomp.blade_mass', ref = 1.e4) - elif opt_options['merit_figure'] == 'LCOE': - wt_opt.model.add_objective('financese.lcoe', ref = 0.1) - elif opt_options['merit_figure'] == 'blade_tip_deflection': - wt_opt.model.add_objective('tcons.tip_deflection_ratio') - elif opt_options['merit_figure'] == 'tower_mass': - wt_opt.model.add_objective('towerse.tower_mass') - elif opt_options['merit_figure'] == 'tower_cost': - wt_opt.model.add_objective('tcc.tower_cost') - elif opt_options['merit_figure'] == 'Cp': - if modeling_options['Analysis_Flags']['ServoSE']: - wt_opt.model.add_objective('sse.powercurve.Cp_regII', ref = -1.) - else: - wt_opt.model.add_objective('ccblade.CP', ref = -1.) - else: - exit('The merit figure ' + opt_options['merit_figure'] + ' is not supported.') - - # Set optimization design variables. - - if blade_opt_options['aero_shape']['twist']['flag']: - indices = range(2, blade_opt_options['aero_shape']['twist']['n_opt']) - wt_opt.model.add_design_var('blade.opt_var.twist_opt_gain', indices = indices, lower=0., upper=1.) - - chord_options = blade_opt_options['aero_shape']['chord'] - if chord_options['flag']: - indices = range(3, chord_options['n_opt'] - 1) - wt_opt.model.add_design_var('blade.opt_var.chord_opt_gain', indices = indices, lower=chord_options['min_gain'], upper=chord_options['max_gain']) - - if blade_opt_options['aero_shape']['af_positions']['flag']: - n_af = modeling_options['blade']['n_af_span'] - indices = range(blade_opt_options['aero_shape']['af_positions']['af_start'],n_af - 1) - af_pos_init = wt_init['components']['blade']['outer_shape_bem']['airfoil_position']['grid'] - lb_af = np.zeros(n_af) - ub_af = np.zeros(n_af) - for i in range(1,indices[0]): - lb_af[i] = ub_af[i] = af_pos_init[i] - for i in indices: - lb_af[i] = 0.5*(af_pos_init[i-1] + af_pos_init[i]) + step_size - ub_af[i] = 0.5*(af_pos_init[i+1] + af_pos_init[i]) - step_size - lb_af[-1] = ub_af[-1] = 1. - wt_opt.model.add_design_var('blade.opt_var.af_position', indices = indices, lower=lb_af[indices], upper=ub_af[indices]) - - spar_cap_ss_options = blade_opt_options['structure']['spar_cap_ss'] - if spar_cap_ss_options['flag']: - indices = range(1,spar_cap_ss_options['n_opt'] - 1) - wt_opt.model.add_design_var('blade.opt_var.spar_cap_ss_opt_gain', indices = indices, lower=spar_cap_ss_options['min_gain'], upper=spar_cap_ss_options['max_gain']) - - # Only add the pressure side design variables if we do set - # `equal_to_suction` as False in the optimization yaml. - spar_cap_ps_options = blade_opt_options['structure']['spar_cap_ps'] - if spar_cap_ps_options['flag'] and not spar_cap_ps_options['equal_to_suction']: - indices = range(1, spar_cap_ps_options['n_opt'] - 1) - wt_opt.model.add_design_var('blade.opt_var.spar_cap_ps_opt_gain', indices = indices, lower=spar_cap_ps_options['min_gain'], upper=spar_cap_ps_options['max_gain']) - - if tower_opt_options['outer_diameter']['flag']: - wt_opt.model.add_design_var('tower.diameter', lower=tower_opt_options['outer_diameter']['lower_bound'], upper=tower_opt_options['outer_diameter']['upper_bound'], ref=5.) - - if tower_opt_options['layer_thickness']['flag']: - wt_opt.model.add_design_var('tower.layer_thickness', lower=tower_opt_options['layer_thickness']['lower_bound'], upper=tower_opt_options['layer_thickness']['upper_bound'], ref=1e-2) - - # -- Control -- - if control_opt_options['tsr']['flag']: - wt_opt.model.add_design_var('opt_var.tsr_opt_gain', lower=control_opt_options['tsr']['min_gain'], - upper=control_opt_options['tsr']['max_gain']) - if control_opt_options['servo']['pitch_control']['flag']: - wt_opt.model.add_design_var('control.PC_omega', lower=control_opt_options['servo']['pitch_control']['omega_min'], - upper=control_opt_options['servo']['pitch_control']['omega_max']) - wt_opt.model.add_design_var('control.PC_zeta', lower=control_opt_options['servo']['pitch_control']['zeta_min'], - upper=control_opt_options['servo']['pitch_control']['zeta_max']) - if control_opt_options['servo']['torque_control']['flag']: - wt_opt.model.add_design_var('control.VS_omega', lower=control_opt_options['servo']['torque_control']['omega_min'], - upper=control_opt_options['servo']['torque_control']['omega_max']) - wt_opt.model.add_design_var('control.VS_zeta', lower=control_opt_options['servo']['torque_control']['zeta_min'], - upper=control_opt_options['servo']['torque_control']['zeta_max']) - - # Set non-linear constraints - blade_constraints = opt_options['constraints']['blade'] - if blade_constraints['strains_spar_cap_ss']['flag']: - if blade_opt_options['structure']['spar_cap_ss']['flag']: - wt_opt.model.add_constraint('rlds.constr.constr_max_strainU_spar', upper= 1.0) - else: - print('WARNING: the strains of the suction-side spar cap are set to be constrained, but spar cap thickness is not an active design variable. The constraint is not enforced.') - - if blade_constraints['strains_spar_cap_ps']['flag']: - if blade_opt_options['structure']['spar_cap_ps']['flag'] or blade_opt_options['structure']['spar_cap_ps']['equal_to_suction']: - wt_opt.model.add_constraint('rlds.constr.constr_max_strainL_spar', upper= 1.0) - else: - print('WARNING: the strains of the pressure-side spar cap are set to be constrained, but spar cap thickness is not an active design variable. The constraint is not enforced.') - - if blade_constraints['stall']['flag']: - if blade_opt_options['aero_shape']['twist']['flag']: - wt_opt.model.add_constraint('stall_check.no_stall_constraint', upper= 1.0) - else: - print('WARNING: the margin to stall is set to be constrained, but twist is not an active design variable. The constraint is not enforced.') - - if blade_constraints['tip_deflection']['flag']: - if blade_opt_options['structure']['spar_cap_ss']['flag'] or blade_opt_options['structure']['spar_cap_ps']['flag']: - wt_opt.model.add_constraint('tcons.tip_deflection_ratio', upper= blade_constraints['tip_deflection']['ratio']) - else: - print('WARNING: the tip deflection is set to be constrained, but spar caps thickness is not an active design variable. The constraint is not enforced.') - - if blade_constraints['chord']['flag']: - if blade_opt_options['aero_shape']['chord']['flag']: - wt_opt.model.add_constraint('blade.pa.max_chord_constr', upper= 1.0) - else: - print('WARNING: the max chord is set to be constrained, but chord is not an active design variable. The constraint is not enforced.') - - if blade_constraints['frequency']['flap_above_3P']: - if blade_opt_options['structure']['spar_cap_ss']['flag'] or blade_opt_options['structure']['spar_cap_ps']['flag']: - wt_opt.model.add_constraint('rlds.constr.constr_flap_f_margin', upper= 0.0) - else: - print('WARNING: the blade flap frequencies are set to be constrained, but spar caps thickness is not an active design variable. The constraint is not enforced.') - - if blade_constraints['frequency']['edge_above_3P']: - wt_opt.model.add_constraint('rlds.constr.constr_edge_f_margin', upper= 0.0) - - if blade_constraints['rail_transport']['flag']: - if blade_constraints['rail_transport']['8_axle']: - wt_opt.model.add_constraint('elastic.rail.constr_LV_8axle_horiz', lower = 0.8, upper= 1.0) - wt_opt.model.add_constraint('elastic.rail.constr_strainPS', upper= 1.0) - wt_opt.model.add_constraint('elastic.rail.constr_strainSS', upper= 1.0) - elif blade_constraints['rail_transport']['4_axle']: - wt_opt.model.add_constraint('elastic.rail.constr_LV_4axle_horiz', upper= 1.0) - else: - exit('You have activated the rail transport constraint module. Please define whether you want to model 4- or 8-axle flatcars.') - - if opt_options['constraints']['blade']['moment_coefficient']['flag']: - wt_opt.model.add_constraint('ccblade.CM', lower= opt_options['constraints']['blade']['moment_coefficient']['min'], upper= opt_options['constraints']['blade']['moment_coefficient']['max']) - if opt_options['constraints']['blade']['match_cl_cd']['flag_cl'] or opt_options['constraints']['blade']['match_cl_cd']['flag_cd']: - data_target = np.loadtxt(opt_options['constraints']['blade']['match_cl_cd']['filename']) - eta_opt = np.linspace(0., 1., opt_options['optimization_variables']['blade']['aero_shape']['twist']['n_opt']) - target_cl = np.interp(eta_opt, data_target[:,0], data_target[:,3]) - target_cd = np.interp(eta_opt, data_target[:,0], data_target[:,4]) - eps_cl = 1.e-2 - if opt_options['constraints']['blade']['match_cl_cd']['flag_cl']: - wt_opt.model.add_constraint('ccblade.cl_n_opt', lower = target_cl-eps_cl, upper = target_cl+eps_cl) - if opt_options['constraints']['blade']['match_cl_cd']['flag_cd']: - wt_opt.model.add_constraint('ccblade.cd_n_opt', lower = target_cd-eps_cl, upper = target_cd+eps_cl) - if opt_options['constraints']['blade']['match_L_D']['flag_L'] or opt_options['constraints']['blade']['match_L_D']['flag_D']: - data_target = np.loadtxt(opt_options['constraints']['blade']['match_L_D']['filename']) - eta_opt = np.linspace(0., 1., opt_options['optimization_variables']['blade']['aero_shape']['twist']['n_opt']) - target_L = np.interp(eta_opt, data_target[:,0], data_target[:,7]) - target_D = np.interp(eta_opt, data_target[:,0], data_target[:,8]) - eps_L = 1.e+2 - if opt_options['constraints']['blade']['match_L_D']['flag_L']: - wt_opt.model.add_constraint('ccblade.L_n_opt', lower = target_L-eps_L, upper = target_L+eps_L) - if opt_options['constraints']['blade']['match_L_D']['flag_D']: - wt_opt.model.add_constraint('ccblade.D_n_opt', lower = target_D-eps_L, upper = target_D+eps_L) - - tower_constraints = opt_options['constraints']['tower'] - if tower_constraints['height_constraint']['flag']: - wt_opt.model.add_constraint('towerse.height_constraint', - lower=tower_constraints['height_constraint']['lower_bound'], - upper=tower_constraints['height_constraint']['upper_bound']) - - if tower_constraints['stress']['flag']: - wt_opt.model.add_constraint('towerse.post.stress', upper=1.0) - - if tower_constraints['global_buckling']['flag']: - wt_opt.model.add_constraint('towerse.post.global_buckling', upper=1.0) - - if tower_constraints['shell_buckling']['flag']: - wt_opt.model.add_constraint('towerse.post.shell_buckling', upper=1.0) - - if tower_constraints['constr_d_to_t']['flag']: - wt_opt.model.add_constraint('towerse.constr_d_to_t', upper=0.0) - - if tower_constraints['constr_taper']['flag']: - wt_opt.model.add_constraint('towerse.constr_taper', lower=0.0) - - if tower_constraints['slope']['flag']: - wt_opt.model.add_constraint('towerse.slope', upper=1.0) - - if tower_constraints['frequency_1']['flag']: - wt_opt.model.add_constraint('towerse.tower.f1', - lower=tower_constraints['frequency_1']['lower_bound'], - upper=tower_constraints['frequency_1']['upper_bound']) - - control_constraints = opt_options['constraints']['control'] - - # Set recorder on the OpenMDAO driver level using the `optimization_log` - # filename supplied in the optimization yaml - if opt_options['recorder']['flag']: - recorder = om.SqliteRecorder(os.path.join(folder_output, opt_options['recorder']['file_name'])) - wt_opt.driver.add_recorder(recorder) - wt_opt.add_recorder(recorder) - - wt_opt.driver.recording_options['excludes'] = ['*_df'] - wt_opt.driver.recording_options['record_constraints'] = True - wt_opt.driver.recording_options['record_desvars'] = True - wt_opt.driver.recording_options['record_objectives'] = True + wt_opt = myopt.set_driver(wt_opt) + wt_opt = myopt.set_objective(wt_opt) + wt_opt = myopt.set_design_variables(wt_opt, wt_init) + wt_opt = myopt.set_constraints(wt_opt) + wt_opt = myopt.set_recorders(wt_opt) # Setup openmdao problem wt_opt.setup() # Load initial wind turbine data from wt_initial to the openmdao problem wt_opt = yaml2openmdao(wt_opt, modeling_options, wt_init) - wt_opt['blade.opt_var.s_opt_twist'] = np.linspace(0., 1., blade_opt_options['aero_shape']['twist']['n_opt']) - if blade_opt_options['aero_shape']['twist']['flag']: - init_twist_opt = np.interp(wt_opt['blade.opt_var.s_opt_twist'], wt_init['components']['blade']['outer_shape_bem']['twist']['grid'], wt_init['components']['blade']['outer_shape_bem']['twist']['values']) - lb_twist = np.array(blade_opt_options['aero_shape']['twist']['lower_bound']) - ub_twist = np.array(blade_opt_options['aero_shape']['twist']['upper_bound']) - wt_opt['blade.opt_var.twist_opt_gain'] = (init_twist_opt - lb_twist) / (ub_twist - lb_twist) - if max(wt_opt['blade.opt_var.twist_opt_gain']) > 1. or min(wt_opt['blade.opt_var.twist_opt_gain']) < 0.: - print('Warning: the initial twist violates the upper or lower bounds of the twist design variables.') - - blade_constraints = opt_options['constraints']['blade'] - wt_opt['blade.opt_var.s_opt_chord'] = np.linspace(0., 1., blade_opt_options['aero_shape']['chord']['n_opt']) - wt_opt['blade.ps.s_opt_spar_cap_ss'] = np.linspace(0., 1., blade_opt_options['structure']['spar_cap_ss']['n_opt']) - wt_opt['blade.ps.s_opt_spar_cap_ps'] = np.linspace(0., 1., blade_opt_options['structure']['spar_cap_ps']['n_opt']) - wt_opt['rlds.constr.max_strainU_spar'] = blade_constraints['strains_spar_cap_ss']['max'] - wt_opt['rlds.constr.max_strainL_spar'] = blade_constraints['strains_spar_cap_ps']['max'] - wt_opt['stall_check.stall_margin'] = blade_constraints['stall']['margin'] * 180. / np.pi + wt_opt = myopt.set_initial(wt_opt, wt_init) # If the user provides values in this dict, they overwrite # whatever values have been set by the yaml files. @@ -370,52 +89,7 @@ def run_wisdem(fname_wt_input, fname_modeling_options, fname_opt_options, overri # Place the last design variables from a previous run into the problem. # This needs to occur after the above setup() and yaml2openmdao() calls # so these values are correctly placed in the problem. - if 'warmstart_file' in opt_options['driver']: - - # Directly read the pyoptsparse sqlite db file - from pyoptsparse import SqliteDict - db = SqliteDict(opt_options['driver']['warmstart_file']) - - # Grab the last iteration's design variables - last_key = db['last'] - desvars = db[last_key]['xuser'] - - # Obtain the already-setup OM problem's design variables - if wt_opt.model._static_mode: - design_vars = wt_opt.model._static_design_vars - else: - design_vars = wt_opt.model._design_vars - - # Get the absolute names from the promoted names within the OM model. - # We need this because the pyoptsparse db has the absolute names for - # variables but the OM model uses the promoted names. - prom2abs = wt_opt.model._var_allprocs_prom2abs_list['output'] - abs2prom = {} - for key in design_vars: - abs2prom[prom2abs[key][0]] = key - - # Loop through each design variable - for key in desvars: - prom_key = abs2prom[key] - - # Scale each DV based on the OM scaling from the problem. - # This assumes we're running the same problem with the same scaling - scaler = design_vars[prom_key]['scaler'] - adder = design_vars[prom_key]['adder'] - - if scaler is None: - scaler = 1.0 - if adder is None: - adder = 0.0 - - scaled_dv = desvars[key] / scaler - adder - - # Special handling for blade twist as we only have the - # last few control points as design variables - if 'twist_opt_gain' in key: - wt_opt[key][2:] = scaled_dv - else: - wt_opt[key][:] = scaled_dv + wt_opt = myopt.set_restart(wt_opt) if 'check_totals' in opt_options['driver']: if opt_options['driver']['check_totals']: diff --git a/wisdem/rotorse/rotor_loads_defl_strains.py b/wisdem/rotorse/rotor_loads_defl_strains.py index ddaa1b604..37ba9dc8f 100644 --- a/wisdem/rotorse/rotor_loads_defl_strains.py +++ b/wisdem/rotorse/rotor_loads_defl_strains.py @@ -316,7 +316,8 @@ def rotate(x,y): lump = 0 # 0= consistent mass matrix, 1= lumped mass matrix tol = 1e-9 # frequency convergence tolerance shift = 0.0 # frequency shift-factor for rigid body modes, make 0 for pos.def. [K] - blade.enableDynamics(self.n_freq, Mmethod, lump, tol, shift) + # Run twice the number of modes to ensure that we can ignore the torsional modes and still get the desired number of fore-aft, side-side modes + blade.enableDynamics(2*self.n_freq, Mmethod, lump, tol, shift) # ---------------------------- # ------ load case 1, blade 1 ------------ @@ -358,8 +359,13 @@ def rotate(x,y): dz = displacements.dz[iCase,:] # Mode shapes and frequencies - freq_x, freq_y, mshapes_x, mshapes_y = util.get_mode_shapes(r, modal.freq, modal.xdsp, modal.ydsp, modal.zdsp, modal.xmpf, modal.ympf, modal.zmpf) - + n_freq2 = int(self.n_freq/2) + freq_x, freq_y, mshapes_x, mshapes_y = util.get_xy_mode_shapes(r, modal.freq, modal.xdsp, modal.ydsp, modal.zdsp, modal.xmpf, modal.ympf, modal.zmpf) + freq_x = freq_x[:n_freq2] + freq_y = freq_y[:n_freq2] + mshapes_x = mshapes_x[:n_freq2,:] + mshapes_y = mshapes_y[:n_freq2,:] + # shear and bending, one per element (convert from local to global c.s.) Fz = np.r_[-forces.Nx[iCase,0], forces.Nx[iCase, 1::2]] Vy = np.r_[-forces.Vy[iCase,0], forces.Vy[iCase, 1::2]] @@ -396,7 +402,7 @@ def strain(xu, yu, xl, yl): # Store outputs outputs['root_F'] = -1.0 * np.array([reactions.Fx.sum(), reactions.Fy.sum(), reactions.Fz.sum()]) outputs['root_M'] = -1.0 * np.array([reactions.Mxx.sum(), reactions.Myy.sum(), reactions.Mzz.sum()]) - outputs['freqs'] = modal.freq + outputs['freqs'] = modal.freq[:self.n_freq] outputs['edge_mode_shapes'] = mshapes_y outputs['flap_mode_shapes'] = mshapes_x # Dense numpy command that interleaves and alternates flap and edge modes diff --git a/wisdem/test/test_commonse/test_utilities.py b/wisdem/test/test_commonse/test_utilities.py index f121dc74d..d404944a8 100644 --- a/wisdem/test/test_commonse/test_utilities.py +++ b/wisdem/test/test_commonse/test_utilities.py @@ -33,6 +33,23 @@ def testModalCoefficients(self): pp = util.get_modal_coefficients(x, y) npt.assert_almost_equal(p[2:]/p[2:].sum(), pp) + def testGetXYModes(self): + r = np.linspace(0,1,20) + n = 10 + n2 = int(n/2) + dx = dy = dz = np.tile(np.r_[r**2 + 10.], (n,1)) + freqs = np.arange(n) + xm = np.zeros(n) + ym = np.zeros(n) + zm = np.zeros(n) + xm[0] = xm[3] = xm[6] = xm[9] = 1 + ym[1] = ym[4] = ym[7] = 1 + zm[2] = zm[5] = zm[8] = 1 + + freq_x, freq_y, _, _ = util.get_xy_mode_shapes(r, freqs, dx, dy, dz, xm, ym, zm) + npt.assert_array_equal(freq_x, np.r_[0, 3, 6, 9, np.zeros(n2-4)]) + npt.assert_array_equal(freq_y, np.r_[1, 4, 7, np.zeros(n2-3)]) + def testRotateI(self): I = np.arange(6)+1 th = np.deg2rad(45)