diff --git a/.gitignore b/.gitignore index f1619701e..4fd190d5d 100644 --- a/.gitignore +++ b/.gitignore @@ -78,6 +78,8 @@ nosetests.xml *.log *.sql *.sqlite +*.err +*.smt # OS generated files # ###################### diff --git a/examples/15_RAFT_Studies/README.md b/examples/15_RAFT_Studies/README.md index 5563d16f7..508dfb1a9 100644 --- a/examples/15_RAFT_Studies/README.md +++ b/examples/15_RAFT_Studies/README.md @@ -22,4 +22,25 @@ To optimize the UMaine semi using raft, run Two options are provided in this example: `../06_IEA-15-240-RWT/IEA-15-240-RWT_VolturnUS-S.yaml` (original) and `raft_opt_out.yaml`, the optimized output. +# RAFT Surrogate Model Training with Design of Experiment +To create the surrogate model with RAFT simulations, run + ``` + python weis_driver_level1_doe.py + ``` + If parallel execution is possible (e.g., HPC) use slurm script to run with MPI across large number of cpus/nodes + ``` + sbatch weis_driver_level1_doe.sh + ``` + or MPI directly for cpus in a local computer + ``` + mpirun -np 4 python weis_driver_level1_doe.py + ``` + + In `analysis_options_level1_doe.yaml`, the design variables (control and plant), driver (optimization:flag: False, design_of_experiments:flag: True), and recorder (flag: True, includes: ['*']) can be set up. Variables marked as design_variables (flag: True) are used as inputs, while all other variables are used as outputs of the surrogate model. Not all parameters are supported at this moment. Seven control parameters, diameters of floating platform members, joint locations, rotor diameter are currently implemented. + + Once design of experiment is completed, a surrogate model for each output variable will be trained (in parallel if MPI is used), and recorder file_name.smt (log_opt.smt in the tutorial case) will be created. + + At this moment, surrogate model training is not stable with large number of samples. + + Further design coupling studies and/or surrogate-based optimization studies can be done in separate script files to be developed. diff --git a/examples/15_RAFT_Studies/analysis_options_level1_doe.yaml b/examples/15_RAFT_Studies/analysis_options_level1_doe.yaml new file mode 100644 index 000000000..18babcb00 --- /dev/null +++ b/examples/15_RAFT_Studies/analysis_options_level1_doe.yaml @@ -0,0 +1,136 @@ +general: + folder_output: outputs/level1_doe + fname_output: refturb_output + +design_variables: + control: + ps_percent: #Peak shaving #defalt: 0.85 + flag: False + lower_bound: 0.75 + upper_bound: 1.0 + servo: + pitch_control: + omega: #rad/s #default: 0.2 + flag: True + min: 0.1 + max: 0.7 + zeta: #rad/s #default: 1.0 + flag: True + min: 0.7 + max: 1.5 + Kp_float: #s #default: -10.0 + flag: True + min: -20.0 #-100.0 + max: 0.0 + ptfm_freq: #rad/s default: 0.2 + flag: True + min: 0.1 #0.00001 + max: 1.5 + torque_control: + omega: #rad/s #default: 0.12 + flag: False + min: 0.1 + max: 0.7 + zeta: #rad/s #default: 0.85 + flag: False + min: 0.7 + max: 1.5 + floating: + members: + flag: False + groups: + - names: [main_column] # Platform sizing: diameter (Center Column) + diameter: + lower_bound: 6.0 #default: 10 + upper_bound: 14.0 + constant: True + - names: [column1, column2, column3] # Platform sizing: diameter (Side Columns) + diameter: + lower_bound: 10.5 #11.0 #10.7 #8.5 # default: 12.5 + upper_bound: 14.5 #14.5 #16.5 + constant: True + - names: [Y_pontoon_upper1, Y_pontoon_upper2, Y_pontoon_upper3] # Platform sizing: Pantoon upper diameter + diameter: + lower_bound: 0.71 # default: 0.91 + upper_bound: 1.11 + constant: True + - names: [Y_pontoon_lower1, Y_pontoon_lower2, Y_pontoon_lower3] # Platform sizing: Pantoon lower diameter + diameter: + lower_bound: 6.6148 #7.6148 # default: 9.6148 + upper_bound: 13.6148 #11.6148 + constant: True + joints: + flag: False + z_coordinate: + - names: [main_keel, col1_keel, col2_keel, col3_keel] #Platform sizing: form -20 (keel) to 15 (freeboard) + lower_bound: -24.0 #default -20 + upper_bound: -16.0 + - names: [main_freeboard, col1_freeboard, col2_freeboard, col3_freeboard] #Platform sizing: form -20 (keel) to 15 (freeboard) + lower_bound: 7.0 #default 15 + upper_bound: 21.0 + r_coordinate: + - names: [col1_keel, col1_freeboard, col2_keel, col2_freeboard, col3_keel, col3_freeboard] #Platform sizing: Spacing + lower_bound: 41.57 #51.57 #48.57 #63.57 #42.57 # default: 51.57 + upper_bound: 61.57 #72.57 #63.57 #60.0 + rotor_diameter: + flag: True + minimum: 235.0 #242.23775645 + maximum: 250.0 + tower: # does not work for design coupling analysis + outer_diameter: + flag: False + lower_bound: 6.0 # 6.5[m] + upper_bound: 7.0 + section_height: + flag: False + lower_bound: 13.0 # 5.0[m] I think in default it has 10 sections: 10 *15=150==hub height + upper_bound: 18.0 + +constraints: + control: + Max_PtfmPitch: + flag: True + max: 6.757845 #6.74303 #deg + Std_PtfmPitch: + flag: True + max: 2.0 #deg + nacelle_acceleration: + flag: True + max: 3.2667 # m/s^2 + tower: + height_constraint: + flag: False + lower_bound: 0.1 + upper_bound: 5.0 + frequency_1: + flag: True + lower_bound: 0.35 + upper_bound: 0.51 + slope: + flag: False + +merit_figure: platform_mass # Merit figure of the optimization problem. The options are 'AEP' - 'LCOE' - 'Cp' - 'blade_mass' - 'blade_tip_deflection' + +driver: + optimization: + flag: False # Flag to enable optimization + solver: COBYLA #LN_COBYLA # Optimization solver. Other options are 'SLSQP' - 'CONMIN' + tol: 1.e-3 # Optimality tolerance + max_iter: 300 # Maximum number of iterations (SLSQP) + design_of_experiments: + flag: True # Flag to enable design of experiments + run_parallel: True # Flag to run using parallel processing + seed: 2 + generator: LatinHypercube #Uniform #LatinHypercube #FullFact # Type of input generator. (Uniform) + num_samples: 50 # number of samples for (Uniform only) + levels: 1 + skip_doe_if_results_exist: True + skip_smt_if_results_exist: True + train_surrogate_model: True + +recorder: + flag: True # Flag to activate OpenMDAO recorder + file_name: log_opt.sql # Name of OpenMDAO recorder + includes: ['*'] #['*raft*','*floating*','*platform*'] + + diff --git a/examples/15_RAFT_Studies/modeling_options_level1.yaml b/examples/15_RAFT_Studies/modeling_options_level1.yaml new file mode 100644 index 000000000..ae27da70b --- /dev/null +++ b/examples/15_RAFT_Studies/modeling_options_level1.yaml @@ -0,0 +1,119 @@ +General: + verbosity: False # When set to True, the code prints to screen many infos + openfast_configuration: + use_exe: True + allow_fails: True + fail_value: 9999 + +WISDEM: + RotorSE: + flag: True + spar_cap_ss: Spar_Cap_SS + spar_cap_ps: Spar_Cap_PS + te_ss: TE_reinforcement_SS + te_ps: TE_reinforcement_PS + TowerSE: + flag: True + DriveSE: + flag: True + FloatingSE: + flag: True + # BOS: + # flag: True + +Level3: # Options for WEIS fidelity level 3 = nonlinear time domain + flag: False + simulation: + DT: 0.01 + CompElast: 1 + CompInflow: 1 + CompAero: 2 + CompServo: 1 + CompHydro: 1 + CompSub: 0 + CompMooring: 3 + CompIce: 0 + OutFileFmt: 3 + linearization: + Linearize: False + ElastoDyn: + FlapDOF1: True + FlapDOF2: True + EdgeDOF: True + TeetDOF: False + DrTrDOF: False + GenDOF: True + YawDOF: False + TwFADOF1 : True + TwFADOF2 : True + TwSSDOF1 : True + TwSSDOF2 : True + PtfmSgDOF: True + PtfmSwDOF: True + PtfmHvDOF: True + PtfmRDOF : True + PtfmPDOF : True + PtfmYDOF : True + HydroDyn: + WvLowCOff: 0.15708 + WvHiCOff: 3.2 + WaveSeed1: 123456789 + AddBQuad1: [9.23e5, 0.0, 0.0, 0.0, -8.92e6, 0.0] + AddBQuad2: [0.0, 9.23e5, 0.0, 8.92e6, 0.0, 0.0] + AddBQuad3: [0.0, 0.0, 2.3e6, 0.0, 0.0, 0.0] + AddBQuad4: [0.0, 8.92e6, 0.0, 1.68e10, 0.0, 0.0] + AddBQuad5: [-8.92e6, 0.0, 0.0, 0.0, 1.68e10, 0.0] + AddBQuad6: [0.0, 0.0, 0.0, 0.0, 0.0, 4.8e10] + PotMod: 1 + # WaveMod: 0 + +Level1: + flag: True + potential_model_override: 0 + trim_ballast: 2 + heave_tol: 1 + save_designs: True + +ROSCO: + flag: True + SD_Mode: 0 + PS_Mode: 1 + ps_percent: 0.85 + F_LPFType: 2 + F_NotchType: 2 + Fl_Mode: 2 + tuning_yaml: ../01_aeroelasticse/OpenFAST_models/IEA-15-240-RWT/IEA-15-240-RWT-UMaineSemi/IEA15MW-UMaineSemi.yaml + zeta_pc: [1] + omega_pc: [0.2] + U_pc: [12] + zeta_vs: 0.85 # Torque controller desired damping ratio [-] + omega_vs: 0.12 + twr_freq: 3.2 + ptfm_freq: 0.2 + Kp_float: -10 + +DLC_driver: + metocean_conditions: + wind_speed: [4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24.] + wave_height_NSS: [0.83, 0.88, 0.94, 1.03, 1.16, 1.34, 1.57, 1.86, 2.22, 2.62, 3.07] + wave_period_NSS: [6.9, 6.96, 7.02, 7.12, 7.25, 7.43, 7.66, 7.94, 8.27, 8.63, 9.01] + wave_height_SSS: [6.3, 8, 8, 8.1, 8.5, 8.5, 9.8, 9.8, 9.8, 9.8, 9.9] + wave_period_SSS: [11.5, 12.7, 12.7, 12.8, 13.1, 13.1, 14.1, 14.1, 14.1, 14.1, 14.1] + wave_height1: 6.98 + wave_period1: 11.7 + wave_height50: 10.68 + wave_period50: 14.2 + DLCs: + - DLC: "1.1" + n_seeds: 1 + # - DLC: "1.3" + # n_seeds: 6 + - DLC: "1.4" + - DLC: "1.5" + - DLC: "1.6" + n_seeds: 1 + - DLC: "6.1" + n_seeds: 1 + # - DLC: "6.3" + # n_seeds: 6 + diff --git a/examples/15_RAFT_Studies/weis_driver_level1_doe.py b/examples/15_RAFT_Studies/weis_driver_level1_doe.py new file mode 100644 index 000000000..d733d89e8 --- /dev/null +++ b/examples/15_RAFT_Studies/weis_driver_level1_doe.py @@ -0,0 +1,25 @@ +import os +import time +import sys + +from weis.glue_code.runWEIS import run_weis +from openmdao.utils.mpi import MPI + +run_dir = os.path.dirname( os.path.realpath(__file__) ) + os.sep +fname_wt_input = os.path.join(run_dir,"..","06_IEA-15-240-RWT", "IEA-15-240-RWT_VolturnUS-S.yaml") +fname_modeling_options = run_dir + "modeling_options_level1.yaml" +fname_analysis_options = run_dir + "analysis_options_level1_doe.yaml" +overridden_values = {} + +tt = time.time() +wt_opt, modeling_options, opt_options = run_weis(fname_wt_input, fname_modeling_options, fname_analysis_options, overridden_values) + +if MPI: + rank = MPI.COMM_WORLD.Get_rank() +else: + rank = 0 +if rank == 0: + print('Run time: %f'%(time.time()-tt)) + sys.stdout.flush() + +print('rank = {:}, exiting'.format(rank)) diff --git a/weis/glue_code/runWEIS.py b/weis/glue_code/runWEIS.py index 92e6d9176..826d480b1 100644 --- a/weis/glue_code/runWEIS.py +++ b/weis/glue_code/runWEIS.py @@ -11,6 +11,7 @@ from weis.control.tmd import assign_TMD_values from weis.aeroelasticse.FileTools import save_yaml from wisdem.inputs.validation import simple_types +from weis.surrogate.WTsurrogate import WindTurbineDOE2SM fd_methods = ['SLSQP','SNOPT', 'LD_MMA'] evolutionary_methods = ['DE', 'NSGA2'] @@ -114,6 +115,7 @@ def run_weis(fname_wt_input, fname_modeling_options, fname_opt_options, geometry else: color_i = 0 rank = 0 + max_cores = 1 # make the folder_output relative to the input, if it's a relative path analysis_input_dir = os.path.dirname(opt_options['fname_input_analysis']) @@ -211,13 +213,42 @@ def run_weis(fname_wt_input, fname_modeling_options, fname_opt_options, geometry checks = wt_opt.check_partials(compact_print=True) sys.stdout.flush() + + # DOE and SMT skip logics + SKIP_DRIVER = False + SKIP_SMT = False + sm_filename = os.path.join(folder_output, os.path.splitext(opt_options['recorder']['file_name'])[0] + '.smt') + sql_filename = os.path.join(folder_output, opt_options['recorder']['file_name']) + if MPI and max_cores>1: + sql_filename += '_{:}'.format(rank) + if opt_options['opt_flag'] and opt_options['driver']['design_of_experiments']['flag']: # if DOE enabled + # Skip DOE Driver if SQL file exists + if opt_options['driver']['design_of_experiments']['skip_doe_if_results_exist']: # if DOE skip flag set + if MPI: + doe_file_exist = MPI.COMM_WORLD.gather(os.path.isfile(sql_filename), root=0) + doe_file_exist = MPI.COMM_WORLD.bcast(doe_file_exist, root=0) + else: + doe_file_exist = [os.path.isfile(sql_filename)] + # If sql files exist for all MPI ranks, and no additional sql file exists outside of MPI size, set skip flag + if all(doe_file_exist) and not os.path.isfile(os.path.join(folder_output, opt_options['recorder']['file_name']) + '_{:}'.format(max_cores)): + SKIP_DRIVER = True + if (not MPI) or (MPI and rank == 0): + print('File {:} exists. Skipping the design of experiments.'.format(sql_filename)) + # Skip SM Training if SMT file exists + if opt_options['driver']['design_of_experiments']['skip_smt_if_results_exist']: # if SMT skip flag set + if os.path.isfile(sm_filename): + SKIP_SMT = True + if (not MPI) or (MPI and rank == 0): + print('File {:} exists. Skipping the design of experiments.'.format(sm_filename)) + # Run openmdao problem if opt_options['opt_flag']: - wt_opt.run_driver() + if not SKIP_DRIVER: + wt_opt.run_driver() else: wt_opt.run_model() - if (not MPI) or (MPI and rank == 0): + if ((not MPI) or (MPI and rank == 0)) and (not SKIP_DRIVER): # Save data coming from openmdao to an output yaml file froot_out = os.path.join(folder_output, opt_options['general']['fname_output']) # Remove the fst_vt key from the dictionary and write out the modeling options @@ -255,6 +286,16 @@ def run_weis(fname_wt_input, fname_modeling_options, fname_opt_options, geometry if MPI and color_i < 1000000: MPI.COMM_WORLD.Barrier() + # If design_of_experiment, recorder flag, train_surrogate_model are all True, + # collect sql files and create smt object + if opt_options['opt_flag'] and (not SKIP_SMT): + if opt_options['driver']['design_of_experiments']['flag'] and opt_options['recorder']['flag']: + if opt_options['driver']['design_of_experiments']['train_surrogate_model']: + WTSM = WindTurbineDOE2SM() + WTSM.read_doe(sql_filename, modeling_options, opt_options) # Parallel reading if MPI + WTSM.train_sm() # Parallel training if MPI + WTSM.write_sm(sm_filename) # Saving will be done in rank=0 + if rank == 0: return wt_opt, modeling_options, opt_options else: diff --git a/weis/inputs/analysis_schema.yaml b/weis/inputs/analysis_schema.yaml index 70917ccd9..be18315e7 100644 --- a/weis/inputs/analysis_schema.yaml +++ b/weis/inputs/analysis_schema.yaml @@ -499,7 +499,28 @@ properties: description: Constrain design to one where OpenFAST simulations don't fail_value default: False - + driver: + type: object + default: {} + properties: + design_of_experiments: + type: object + description: Specification of the design of experiments driver parameters + default: {} + properties: + skip_doe_if_results_exist: + type: boolean + default: True + description: Toggle skipping the design of experiments if recorded DOE results (.sql) exist + skip_smt_if_results_exist: + type: boolean + default: True + description: Toggle skipping the surrogate model training if surrogate model file (.smt) exist + train_surrogate_model: + type: boolean + default: False + description: Train surrogate model after running design of experiment + merit_figure: type: string diff --git a/weis/surrogate/WTsurrogate.py b/weis/surrogate/WTsurrogate.py new file mode 100644 index 000000000..53a513c5f --- /dev/null +++ b/weis/surrogate/WTsurrogate.py @@ -0,0 +1,677 @@ +import numpy as np +import csv +import os +import time +import re +import pickle as pkl +import openmdao.api as om +from openmdao.utils.mpi import MPI +from smt.surrogate_models import SGP + +class WindTurbineDOE2SM(): + # Read DOE sql files, Process data, Create and train surrogate models, and Save them in smt file. + + def __init__(self): + self._doe_loaded = False + self._sm_trained = False + + def read_doe(self, sql_file, modeling_options, opt_options): + # Read DOE sql files. If MPI, read them in parallel. + + if MPI: + rank = MPI.COMM_WORLD.Get_rank() + else: + rank = 0 + + cr = om.CaseReader(sql_file) + print('reading cases from {:}'.format(sql_file)) + cases = cr.list_cases('driver') + + if (not MPI) or (MPI and rank == 0): + case = cases[0] + inputs = cr.get_case(case).inputs + outputs = cr.get_case(case).outputs + input_keys_ref = list(set(inputs.keys()).intersection([ + 'floating.member_main_column:outer_diameter', + 'floating.member_column1:outer_diameter', + 'floating.member_column2:outer_diameter', + 'floating.member_column3:outer_diameter', + 'floating.member_Y_pontoon_upper1:outer_diameter', + 'floating.member_Y_pontoon_upper2:outer_diameter', + 'floating.member_Y_pontoon_upper3:outer_diameter', + 'floating.member_Y_pontoon_lower1:outer_diameter', + 'floating.member_Y_pontoon_lower2:outer_diameter', + 'floating.member_Y_pontoon_lower3:outer_diameter', + ])) + input_keys_ref.sort() + input_keys_dv = self._identify_dv(input_keys_ref, opt_options, inputs, outputs) # identify which parameter included in the input_keys_ref are flagged as DV + output_keys_ref = list(set(outputs.keys()).intersection([ + 'tune_rosco_ivc.ps_percent', + 'tune_rosco_ivc.omega_pc', + 'tune_rosco_ivc.zeta_pc', + 'tune_rosco_ivc.Kp_float', + 'tune_rosco_ivc.ptfm_freq', + 'tune_rosco_ivc.omega_vs', + 'tune_rosco_ivc.zeta_vs', + 'configuration.rotor_diameter_user', + 'towerse.tower.fore_aft_freqs', # 123 + 'towerse.tower.side_side_freqs', # 123 + 'towerse.tower.torsion_freqs', # 123 + 'towerse.tower.top_deflection', + 'floating.member_main_column:joint1', # xyz + 'floating.member_main_column:joint2', # xyz + 'floating.member_column1:joint1', # xyz + 'floating.member_column1:joint2', # xyz + 'floating.member_column2:joint1', # xyz + 'floating.member_column2:joint2', # xyz + 'floating.member_column3:joint1', # xyz + 'floating.member_column3:joint2', # xyz + 'floating.jointdv_0', # keel z-location + 'floating.jointdv_1', # freeboard z-location + 'floating.jointdv_2', # column123 r-location + 'raft.Max_Offset', # Maximum distance in surge/sway direction [m] + 'raft.heave_avg', # Average heave over all cases [m] + 'raft.Max_PtfmPitch', # Maximum platform pitch over all cases [deg] + 'raft.Std_PtfmPitch', # Average platform pitch std. over all cases [deg] + 'rigid_body_periods', # Rigid body natural period [s] + 'raft.heave_period', # Heave natural period [s] + 'raft.pitch_period', # Pitch natural period [s] + 'raft.roll_period', # Roll natural period [s] + 'raft.surge_period', # Surge natural period [s] + 'raft.sway_period', # Sway natural period [s] + 'raft.yaw_period', # Yaw natural period [s] + 'raft.max_nac_accel', # Maximum nacelle accelleration over all cases [m/s**2] + 'raft.max_tower_base', # Maximum tower base moment over all cases [N*m] + 'raft.platform_total_center_of_mass', # xyz + 'raft.platform_displacement', + 'raft.platform_mass', # Platform mass + 'tcons.tip_deflection_ratio', # Blade tip deflection ratio (constrained to be <=1.0) + 'financese.lcoe', # WEIS LCOE from FinanceSE + 'rotorse.rp.AEP', # WISDEM AEP from RotorSE + 'rotorse.blade_mass', # Blade mass + #'towerse.tower_mass', # Tower mass + 'fixedse.structural_mass', # System structural mass for fixed foundation turbines + 'floatingse.system_structural_mass', # System structural mass for floating turbines + 'floatingse.platform_mass', # Platform mass from FloatingSE + 'floatingse.platform_cost', # Platform cost + #'floatingse.mooring_mass', # Mooring mass + #'floatingse.mooring_cost', # Mooring cost + 'floatingse.structural_frequencies', + ])) + output_keys_ref.sort() + output_keys_dv = self._identify_dv(output_keys_ref, opt_options, inputs, outputs) # identify which parameter included in the output_keys_ref are flagged as DV + else: + input_keys_ref = None + input_keys_dv = None + output_keys_ref = None + output_keys_dv = None + if MPI: + # broadcast lists of openmdao variable names to all MPI cores + input_keys_ref = MPI.COMM_WORLD.bcast(input_keys_ref, root=0) + input_keys_dv = MPI.COMM_WORLD.bcast(input_keys_dv, root=0) + output_keys_ref = MPI.COMM_WORLD.bcast(output_keys_ref, root=0) + output_keys_dv = MPI.COMM_WORLD.bcast(output_keys_dv, root=0) + + # Retrieve values of the parameters included in the key lists + # from each sql file associated with each MPI core. + dataset = [] + for case in cases: + inputs = cr.get_case(case).inputs + outputs = cr.get_case(case).outputs + var_keys = [] + var_dv = [] + var_values = [] + for key in input_keys_ref: # for each parameter key included in the openmdao inputs object + if len(inputs[key]) == 1: # if scalar + var_keys.append(key) + try: + dvidx = input_keys_dv[input_keys_ref.index(key)] + if ((type(dvidx) == bool) and (dvidx == False)) or \ + ((type(dvidx) == list) and (len(dvidx) == 0)): # If not DV (dvidx is either False or []) + var_dv.append(False) + elif ((type(dvidx) == list) and (len(dvidx) == 1)): # If DV (dvidx is [0], meaning that inputs[key][0] is DV) + var_dv.append(True) + else: + raise Exception + except: + var_dv.append(False) + var_values.append(inputs[key][0]) + else: # if vector + for idx in range(len(np.squeeze(inputs[key]))): + var_keys.append(key + '_' + str(idx)) + try: + dvidx = input_keys_dv[input_keys_ref.index(key)] + if ((type(dvidx) == bool) and (dvidx == False)) or \ + ((type(dvidx) == list) and (len(dvidx) == 0)): # If not DV (dvidx is either False or []) + var_dv.append(False) + elif ((type(dvidx) == list) and (len(dvidx) > 0)): # If DV (dvidx is [i1, i2, ...], meaning that inputs[key][i1], inputs[key][i2], ... are DVs) + for jdx in range(len(dvidx)): + var_dv_append = False + if idx == dvidx[jdx]: + var_dv_append = True + var_dv.append(var_dv_append) # If any element in the vector is DV, then add into the list of DVs + else: + raise Exception + except: + var_dv.append(False) + var_values.append(np.squeeze(inputs[key])[idx]) + for key in output_keys_ref: # for each parameter key included in the openmdao outputs object + if len(outputs[key]) == 1: # if scalar + var_keys.append(key) + try: + dvidx = output_keys_dv[output_keys_ref.index(key)] + if ((type(dvidx) == bool) and (dvidx == False)) or \ + ((type(dvidx) == list) and (len(dvidx) == 0)): # If not DV (dvidx is either False or []) + var_dv.append(False) + elif ((type(dvidx) == list) and (len(dvidx) == 1)): # If DV (dvidx is [0], meaning that outputss[key][0] is DV) + var_dv.append(True) + else: + raise Exception + except: + var_dv.append(False) + var_values.append(outputs[key][0]) + else: # if vector + for idx in range(len(np.squeeze(outputs[key]))): + var_keys.append(key + '_' + str(idx)) + try: + dvidx = output_keys_dv[output_keys_ref.index(key)] + if ((type(dvidx) == bool) and (dvidx == False)) or \ + ((type(dvidx) == list) and (len(dvidx) == 0)): # If not DV (dvidx is either False or []) + var_dv.append(False) + elif ((type(dvidx) == list) and (len(dvidx) > 0)): # If DV (dvidx is [i1, i2, ...], meaning that outputs[key][i1], outputss[key][i2], ... are DVs) + for jdx in range(len(dvidx)): + var_dv_append = False + if idx == dvidx[jdx]: + var_dv_append = True + var_dv.append(var_dv_append) # If any element in the vector is DV, then add into the list of DVs + else: + raise Exception + except: + var_dv.append(False) + var_values.append(np.squeeze(outputs[key])[idx]) + dataset.append(var_values) + + # Gather data values from all MPI cores if running parallel + if MPI: + dataset = MPI.COMM_WORLD.gather(dataset, root=0) + else: + dataset = [dataset] + if rank == 0: + dataset = np.array([dp for proc in dataset for dp in proc]) + # Remove duplicated columns + flag = dataset.shape[1]*[True] + for idx in range(1,dataset.shape[1]): + for jdx in range(idx): + if (np.array_equal(dataset[:,jdx], dataset[:,idx])): + flag[idx] = False + if flag[jdx] == True: + var_keys[jdx] += ('+' + var_keys[idx]) + data_vals = np.zeros(shape=(dataset.shape[0],0)) + data_dv = [] + data_keys = [] + for idx in range(dataset.shape[1]): + if flag[idx]: + data_vals = np.concatenate( + (data_vals, dataset[:,idx].reshape(len(dataset[:,idx]),1)), + axis=1) + data_dv.append(var_dv[idx]) + data_keys.append(var_keys[idx]) + else: + dataset = None + data_vals = None + data_dv = None + data_keys = None + + # Store data + self.data_vals = data_vals + self.data_dv = data_dv + self.data_keys = data_keys + self._doe_loaded = True + + if rank==0: + fname = os.path.join( + opt_options['general']['folder_output'], + os.path.splitext(opt_options['recorder']['file_name'])[0]+'_doedata.pkl' + ) + print('Saving DOE data into: {:}'.format(fname)) + with open(fname, 'wb') as fid: + pkl.dump( + { + 'data_vals': data_vals, + 'data_dv': data_dv, + 'data_keys': data_keys, + }, + fid, protocol=5, + ) + fname = os.path.join( + opt_options['general']['folder_output'], + os.path.splitext(opt_options['recorder']['file_name'])[0]+'_doedata.csv' + ) + print('Saving DOE data into {:}'.format(fname)) + with open(fname, 'wt', newline='') as fid: + dwriter = csv.writer(fid, delimiter=',', quotechar='"', quoting=csv.QUOTE_MINIMAL) + dwriter.writerow(data_keys) + dwriter.writerow(data_dv) + for idx in range(data_vals.shape[0]): + if np.iscomplex(data_vals[idx,:]).any(): + raise Exception('Complex number detected from idx = {:}, key = {:}'.format(idx, data_keys[idx])) + else: + dwriter.writerow(data_vals[idx,:]) + + + + def _identify_dv(self, keys, opt_options, inputs, outputs): + # Identify design variables to set them as input parameters in the surrogate model + # if they are flagged as DV in the analysis_options.yaml file + + dvflag = len(keys)*[False] + + if opt_options['design_variables']['control']['ps_percent']['flag']: + for key in ['tune_rosco_ivc.ps_percent']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['control']['servo']['pitch_control']['omega']['flag']: + for key in ['tune_rosco_ivc.omega_pc']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['control']['servo']['pitch_control']['zeta']['flag']: + for key in ['tune_rosco_ivc.zeta_pc']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['control']['servo']['pitch_control']['Kp_float']['flag']: + for key in ['tune_rosco_ivc.Kp_float']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['control']['servo']['pitch_control']['ptfm_freq']['flag']: + for key in ['tune_rosco_ivc.ptfm_freq']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['control']['servo']['torque_control']['omega']['flag']: + for key in ['tune_rosco_ivc.omega_vs']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['control']['servo']['torque_control']['zeta']['flag']: + for key in ['tune_rosco_ivc.zeta_vs']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['floating']['members']['flag']: + groups = opt_options['design_variables']['floating']['members']['groups'] + for group in groups: + names = group['names'] + for key in keys: + for name in names: + txt = name + ':outer_diameter' + ltxt = len(txt) + if key[-ltxt:] == txt: + try: + idx = keys.index(key) + if group['diameter']['constant']: + jdx = 0 + else: + jdx = True + if not dvflag[idx]: + dvflag[idx] = [jdx] + else: + dvflag[idx].append(jdx) + except: + pass + + if opt_options['design_variables']['floating']['joints']['flag']: + key_prefix = 'floating.jointdv_' + for key in keys: + if key[:len(key_prefix)] == key_prefix: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + if opt_options['design_variables']['rotor_diameter']['flag']: + for key in ['configuration.rotor_diameter_user']: + try: + idx = keys.index(key) + if not dvflag[idx]: + dvflag[idx] = [0] + else: + dvflag[idx].append(0) + except: + pass + + return dvflag + + + def write_sm(self, sm_filename): + # Write trained surrogate models in the smt file in pickle format + + if MPI: + rank = MPI.COMM_WORLD.Get_rank() + n_cores = MPI.COMM_WORLD.Get_size() + else: + rank = 0 + n_cores = 1 + + # File writing only on rank=0 + if rank == 0: + if not self._doe_loaded: + raise Exception('DOE data needs to be loaded first.') + if not self._sm_trained: + raise Exception('SM needs to be trained before saving to file.') + + try: + with open(sm_filename, 'wb') as fid: + pkl.dump(self.dataset_list, fid, protocol=5) + print('Surrogate model written in {:}'.format(sm_filename)) + except: + print('Unable to write surrogate model file: {:}.'.format(sm_filename)) + raise Exception('Unable to write surrogate model file: {:}.'.format(sm_filename)) + + + def train_sm(self): + # Surrogate model training + + if MPI: + rank = MPI.COMM_WORLD.Get_rank() + n_cores = MPI.COMM_WORLD.Get_size() + else: + rank = 0 + n_cores = 1 + + # Prepare dataset columns individually (and split if MPI) + if (not MPI) or (MPI and rank == 0): + data_vals = self.data_vals + data_dv = self.data_dv + data_keys = self.data_keys + + # Number of design variables (n_dv), output variables (n_out), and overall (n_var) + n_dv = data_dv.count(True) + n_out = len(data_dv) - n_dv + n_var = n_dv + n_out + n_data = data_vals.shape[0] + # Lists of indexes for design variables (i_dv) and output variables (i_out) + i_dv = [i for i, e in enumerate(data_dv) if e == True] + i_out = [i for i, e in enumerate(data_dv) if not e == True] + + # Data arrays and their list + data_arr_dv = np.zeros(shape=(n_data, n_dv), dtype=float) + data_lst_out = list(n_out*[None]) + data_dv_keys = list(n_dv*[None]) + data_out_keys = list(n_out*[None]) + data_dv_bounds = list(n_dv*[None]) + data_out_bounds = list(n_out*[None]) + # DVs + for idx in range(len(i_dv)): + jdx = i_dv[idx] + vals = data_vals[:,jdx] + bounds = [min(vals), max(vals)] + if ((bounds[1] - bounds[0]) < 10.0*np.finfo(np.float64).eps): + bounds[1] = 1.0 + bounds[0] = 0.0 + data_arr_dv[:,idx] = (vals - bounds[0])/(bounds[1] - bounds[0]) + data_dv_keys[idx] = data_keys[jdx] + data_dv_bounds[idx] = bounds + data_dv_bounds = np.array(data_dv_bounds).T + # Output variables + for idx in range(len(i_out)): + jdx = i_out[idx] + vals = np.array(data_vals[:,jdx]).reshape(n_data, 1) + bounds = [min(vals)[0], max(vals)[0]] + if ((bounds[1] - bounds[0]) < 10.0*np.finfo(np.float64).eps): + bounds[1] = 1.0 + bounds[0] = 0.0 + data_lst_out[idx] = (vals - bounds[0])/(bounds[1] - bounds[0]) + data_out_keys[idx] = data_keys[jdx] + data_out_bounds[idx] = np.array(bounds) + + # Construct data structure for parallel distribution + dataset_list = [] + for idx in range(len(data_out_keys)): + data_entry = { + 'inputs': { + 'keys': data_dv_keys, + 'vals': data_arr_dv, + 'bounds': data_dv_bounds, + }, + 'outputs': { + 'keys': data_out_keys[idx], + 'vals': data_lst_out[idx], + 'bounds': data_out_bounds[idx], + }, + 'surrogate': None + } + dataset_list.append(data_entry) + + # Distribute for parallel training + if MPI: + dataset_lists = list(self._split_list_chunks(dataset_list, n_cores)) + if len(dataset_lists) < n_cores: + for _ in range(0, n_cores - len(dataset_lists)): + dataset_lists.append([]) + else: + dataset_lists = [dataset_list] + else: + dataset_list = [] + dataset_lists = [] + # Scatter data to train + if MPI: + MPI.COMM_WORLD.barrier() + dataset_list = MPI.COMM_WORLD.scatter(dataset_lists, root=0) + + # Train SM + for data_entry in dataset_list: + t = time.time() + try: + rng = np.random.RandomState(42) + outvalarr = data_entry['outputs']['vals'] + outvalavg = np.mean(outvalarr) + # If response is constant + if np.all(np.isclose(outvalarr, outvalavg, atol=np.abs(outvalavg)*1.0e-3)): + outvalconst = True + R_squared = 1.0 + else: + outvalconst = False + n_data = data_entry['inputs']['vals'].shape[0] + # If number of sample data is enough + if n_data > 30: + n_data_80 = int(0.8*float(n_data)) + # Train the validation SM + smval = SGP_WT(eval_noise=True, print_global=False) + smval.set_training_values( + data_entry['inputs']['vals'][:n_data_80,:], data_entry['outputs']['vals'][:n_data_80,:]) + n_inducing = int(np.min([n_data_80/1.5, 100])) + random_idx = rng.permutation(n_data_80)[:n_inducing] + Z = data_entry['inputs']['vals'][random_idx].copy() + smval.set_inducing_inputs(Z=Z) + smval.train() + smval.set_bounds( + data_entry['inputs']['bounds'], data_entry['outputs']['bounds']) + smval.keys_in = data_entry['inputs']['keys'] + for idx in range(len(smval.keys_in)): + s = smval.keys_in[idx] + i = re.search('[^+]*', s).span() + smval.keys_in[idx] = s[i[0]:i[1]] + # Validate SM + x_val = data_entry['inputs']['vals'][n_data_80:,:] + y_val_true = data_entry['outputs']['vals'][n_data_80:,:] + lb_val_in = np.tile(smval.bounds_in[0,:], (x_val.shape[0], 1)) + ub_val_in = np.tile(smval.bounds_in[1,:], (x_val.shape[0], 1)) + lb_val_out = np.tile(smval.bounds_out[0], (x_val.shape[0], 1)) + ub_val_out = np.tile(smval.bounds_out[1], (x_val.shape[0], 1)) + x_val = lb_val_in + (ub_val_in - lb_val_in)*x_val + y_val_true = lb_val_out + (ub_val_out - lb_val_out)*y_val_true + y_val_pred, v_val_pred = smval.predict(x_val) + SSE = np.sum((y_val_pred - y_val_true)**2) + SST = np.sum((y_val_true - np.mean(y_val_true))**2) + R_squared = 1.0 - SSE/SST + print('rank {:}, R_squared: {:}'.format(rank, R_squared)) + if R_squared <= 0.25: + outvalconst = True + outvalavg = np.mean(outvalarr) + R_squared = 0.0 + else: + R_squared = None + # Train the full SM + data_entry['surrogate'] = SGP_WT(eval_noise=True, print_global=False) + data_entry['surrogate'].set_training_values( + data_entry['inputs']['vals'], data_entry['outputs']['vals']) + n_inducing = int(np.min([n_data/1.5, 100])) + random_idx = rng.permutation(n_data)[:n_inducing] + Z = data_entry['inputs']['vals'][random_idx].copy() + if outvalconst: + data_entry['surrogate'].constant = True + data_entry['surrogate'].constant_value = outvalavg + else: + data_entry['surrogate'].set_inducing_inputs(Z=Z) + data_entry['surrogate'].train() + data_entry['surrogate'].set_bounds( + data_entry['inputs']['bounds'], data_entry['outputs']['bounds']) + data_entry['surrogate'].keys_in = data_entry['inputs']['keys'] + for idx in range(len(data_entry['surrogate'].keys_in)): + s = data_entry['surrogate'].keys_in[idx] + i = re.search('[^+]*', s).span() + data_entry['surrogate'].keys_in[idx] = s[i[0]:i[1]] + data_entry['surrogate'].keys_out = data_entry['outputs']['keys'] + s = data_entry['surrogate'].keys_out + i = re.search('[^+]*', s).span() + data_entry['surrogate'].keys_out = s[i[0]:i[1]] + data_entry['surrogate'].R_squared = R_squared + except: + print('rank {:}, Surrogate model training failed.'.format(rank)) + raise Exception('rank {:}, Surrogate model training failed.'.format(rank)) + t = time.time() - t + print('rank {:}, Surrogate model training done. Time (sec): {:}'.format(rank, t)) + + + # Gather trained data + if MPI: + MPI.COMM_WORLD.barrier() + dataset_lists = MPI.COMM_WORLD.gather(dataset_list, root=0) + if rank == 0: + dataset_list = [item for items in dataset_lists for item in items] + else: + dataset_list = [] + + self.dataset_list = dataset_list + self._sm_trained = True + + + + def _split_list_chunks(self, fulllist, max_n_chunk=1, item_count=None): + # Split items in a list into nested multiple (max_n_chunk) lists in an outer list + # This method is useful to divide jobs for parallel surrogate model training as the number + # of surrogate models trained are generally not equal to the number of parallel cores used. + + item_count = item_count or len(fulllist) + n_chunks = min(item_count, max_n_chunk) + fulllist = iter(fulllist) + floor = item_count // n_chunks + ceiling = floor + 1 + stepdown = item_count % n_chunks + for x_i in range(n_chunks): + length = ceiling if x_i < stepdown else floor + yield [next(fulllist) for _ in range(length)] + + +class SGP_WT(SGP): + # Sparse Gaussian Process surrogate model class, specifically tailored for WEIS use, + # inherited from SGP class included in the SMT Toolbox. + + def _initialize(self): + super()._initialize() + self._bounds_set = False + self.keys_in = None + self.keys_out = None + self.R_squared = None + self.constant = False + self.constant_value = 0.0 + + def set_bounds(self, bounds_in, bounds_out): + self.bounds_in = bounds_in + self.bounds_out = bounds_out + self._bounds_set = True + + def predict(self, x): + # Predicts surrogate model response and variance + # Input (x) and outputs (y_out, v_out) are denormalized (raw scale) values + # while actual surrogate model computation is done in normalized scale. + + x_in = np.array(x) + + if not self._bounds_set: + raise Exception('Normalizing bounds are needed before accessing surrogate model.') + + if not len(x_in.shape) == 2: + raise Exception('Input array x needs to have shape = (:,n_dv).') + + lb_in = np.tile(self.bounds_in[0,:], (x_in.shape[0], 1)) + ub_in = np.tile(self.bounds_in[1,:], (x_in.shape[0], 1)) + + lb_out = np.tile(self.bounds_out[0], (x_in.shape[0], 1)) + ub_out = np.tile(self.bounds_out[1], (x_in.shape[0], 1)) + + x_in_normalized = (x_in - lb_in)/(ub_in - lb_in) + if self.constant: + y_out_normalized = self.constant_value*np.ones((x_in.shape[0], 1), dtype=float) + sqrt_v_out_normalized = 0.0 + else: + y_out_normalized = self.predict_values(x_in_normalized) + sqrt_v_out_normalized = np.sqrt(self.predict_variances(x_in_normalized)) + + y_out = lb_out + (ub_out - lb_out)*y_out_normalized + sqrt_v_out = (ub_out - lb_out)*sqrt_v_out_normalized + v_out = sqrt_v_out**2 + + return y_out, v_out # values and variances + + diff --git a/weis/surrogate/__init__.py b/weis/surrogate/__init__.py new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/weis/surrogate/__init__.py @@ -0,0 +1 @@ +