diff --git a/.gitignore b/.gitignore index 0b2f951ce..e274d706c 100644 --- a/.gitignore +++ b/.gitignore @@ -109,6 +109,7 @@ venv/ ENV/ env.bak/ venv.bak/ +.github/*_environment.yml # Spyder project settings .spyderproject diff --git a/aviary/interface/default_phase_info/two_dof_fiti.py b/aviary/interface/default_phase_info/two_dof_fiti.py index 75f2e1a25..29086e56d 100644 --- a/aviary/interface/default_phase_info/two_dof_fiti.py +++ b/aviary/interface/default_phase_info/two_dof_fiti.py @@ -54,6 +54,7 @@ def create_2dof_based_ascent_phases( ) accel_vals = {} + # need to set trigger units climb1_kwargs = dict( input_speed_type=SpeedType.EAS, input_speed_units='kn', diff --git a/aviary/interface/methods_for_level2.py b/aviary/interface/methods_for_level2.py index 38c2866ac..601dc045d 100644 --- a/aviary/interface/methods_for_level2.py +++ b/aviary/interface/methods_for_level2.py @@ -302,7 +302,7 @@ def load_inputs(self, aviary_inputs, phase_info=None, engine_builder=None): elif self.mission_method is SOLVED: from aviary.interface.default_phase_info.solved import phase_info - print('Loaded default phase_info for' + print('Loaded default phase_info for ' f'{self.mission_method.value.lower()} equations of motion') # create a new dictionary that only contains the phases from phase_info diff --git a/aviary/mission/flops_based/ode/landing_ode.py b/aviary/mission/flops_based/ode/landing_ode.py index c5e6d293d..4681acf1b 100644 --- a/aviary/mission/flops_based/ode/landing_ode.py +++ b/aviary/mission/flops_based/ode/landing_ode.py @@ -12,10 +12,12 @@ from aviary.mission.flops_based.ode.landing_eom import FlareEOM, StallSpeed from aviary.mission.flops_based.ode.takeoff_ode import TakeoffODE as _TakeoffODE from aviary.mission.gasp_based.flight_conditions import FlightConditions +from aviary.mission.gasp_based.ode.time_integration_base_classes import add_SGM_required_inputs from aviary.utils.aviary_values import AviaryValues from aviary.utils.functions import set_aviary_initial_values, promote_aircraft_and_mission_vars from aviary.variable_info.variables import Aircraft, Dynamic, Mission from aviary.variable_info.variables_in import VariablesIn +from aviary.variable_info.enums import AnalysisScheme class ExternalSubsystemGroup(om.Group): @@ -60,14 +62,29 @@ def initialize(self): 'external_subsystems', default=[], desc='list of external subsystem builder instances to be added to the ODE') + self.options.declare( + "analysis_scheme", + default=AnalysisScheme.COLLOCATION, + types=AnalysisScheme, + desc="The analysis method that will be used to close the trajectory; for example collocation or time integration", + ) + def setup(self): options = self.options nn = options["num_nodes"] + analysis_scheme = options['analysis_scheme'] aviary_options = options['aviary_options'] subsystem_options = options['subsystem_options'] core_subsystems = options['core_subsystems'] + if analysis_scheme is AnalysisScheme.SHOOTING: + SGM_required_inputs = { + 't_curr': {'units': 's'}, + Dynamic.Mission.DISTANCE: {'units': 'm'}, + } + add_SGM_required_inputs(self, SGM_required_inputs) + self.add_subsystem( 'input_port', VariablesIn(aviary_options=aviary_options), promotes_inputs=['*'], promotes_outputs=['*']) diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index 71fe5b932..87756a921 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -3,11 +3,13 @@ from dymos.models.atmosphere import USatm1976Comp from aviary.mission.flops_based.ode.mission_EOM import MissionEOM +from aviary.mission.gasp_based.ode.time_integration_base_classes import add_SGM_required_inputs, add_SGM_required_outputs from aviary.utils.aviary_values import AviaryValues from aviary.utils.functions import promote_aircraft_and_mission_vars from aviary.variable_info.variable_meta_data import _MetaData -from aviary.variable_info.variables import Dynamic, Mission +from aviary.variable_info.variables import Dynamic, Mission, Aircraft from aviary.variable_info.variables_in import VariablesIn +from aviary.variable_info.enums import AnalysisScheme class ExternalSubsystemGroup(om.Group): @@ -44,15 +46,29 @@ def initialize(self): values=['path_constraint', 'boundary_constraint', 'bounded', None], desc='flag to enforce throttle constraints on the path or at the segment boundaries or using solver bounds' ) + self.options.declare( + "analysis_scheme", + default=AnalysisScheme.COLLOCATION, + types=AnalysisScheme, + desc="The analysis method that will be used to close the trajectory; for example collocation or time integration", + ) def setup(self): options = self.options nn = options['num_nodes'] + analysis_scheme = options['analysis_scheme'] aviary_options = options['aviary_options'] core_subsystems = options['core_subsystems'] subsystem_options = options['subsystem_options'] engine_count = len(aviary_options.get_val('engine_models')) + if analysis_scheme is AnalysisScheme.SHOOTING: + SGM_required_inputs = { + 't_curr': {'units': 's'}, + Dynamic.Mission.DISTANCE: {'units': 'm'}, + } + add_SGM_required_inputs(self, SGM_required_inputs) + self.add_subsystem( 'input_port', VariablesIn(aviary_options=aviary_options, @@ -169,9 +185,12 @@ def setup(self): promotes_inputs=['*'], promotes_outputs=['*']) + self.set_input_defaults(Dynamic.Mission.MACH, val=np.ones(nn), units='unitless') self.set_input_defaults(Dynamic.Mission.MASS, val=np.ones(nn), units='kg') self.set_input_defaults(Dynamic.Mission.VELOCITY, val=np.ones(nn), units='m/s') self.set_input_defaults(Dynamic.Mission.ALTITUDE, val=np.ones(nn), units='m') + self.set_input_defaults(Dynamic.Mission.ALTITUDE_RATE, + val=np.ones(nn), units='m/s') self.set_input_defaults( Dynamic.Mission.THROTTLE, val=np.ones((nn, engine_count)), units='unitless') @@ -198,6 +217,14 @@ def setup(self): ], promotes_outputs=['initial_mass_residual']) + if analysis_scheme is AnalysisScheme.SHOOTING: + SGM_required_outputs = { + Dynamic.Mission.ALTITUDE_RATE: {'units': 'm/s'}, + } + add_SGM_required_outputs(self, SGM_required_outputs) + + print_level = 0 if analysis_scheme is AnalysisScheme.SHOOTING else 2 + self.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, atol=1.0e-10, rtol=1.0e-10, @@ -205,4 +232,4 @@ def setup(self): self.nonlinear_solver.linesearch = om.BoundsEnforceLS() self.linear_solver = om.DirectSolver(assemble_jac=True) self.nonlinear_solver.options['err_on_non_converge'] = True - self.nonlinear_solver.options['iprint'] = 2 + self.nonlinear_solver.options['iprint'] = print_level diff --git a/aviary/mission/flops_based/ode/takeoff_ode.py b/aviary/mission/flops_based/ode/takeoff_ode.py index 23a0325bb..0c7396a00 100644 --- a/aviary/mission/flops_based/ode/takeoff_ode.py +++ b/aviary/mission/flops_based/ode/takeoff_ode.py @@ -7,10 +7,12 @@ from aviary.mission.flops_based.ode.takeoff_eom import StallSpeed, TakeoffEOM from aviary.mission.gasp_based.flight_conditions import FlightConditions +from aviary.mission.gasp_based.ode.time_integration_base_classes import add_SGM_required_inputs from aviary.utils.aviary_values import AviaryValues from aviary.utils.functions import set_aviary_initial_values, promote_aircraft_and_mission_vars from aviary.variable_info.variables import Aircraft, Dynamic, Mission from aviary.variable_info.variables_in import VariablesIn +from aviary.variable_info.enums import AnalysisScheme class ExternalSubsystemGroup(om.Group): @@ -58,14 +60,29 @@ def initialize(self): 'climbing', default=False, types=bool, desc='mode of operation (ground roll or flight)') + self.options.declare( + "analysis_scheme", + default=AnalysisScheme.COLLOCATION, + types=AnalysisScheme, + desc="The analysis method that will be used to close the trajectory; for example collocation or time integration", + ) + def setup(self): options = self.options nn = options["num_nodes"] + analysis_scheme = options['analysis_scheme'] aviary_options = options['aviary_options'] subsystem_options = options['subsystem_options'] core_subsystems = options['core_subsystems'] + if analysis_scheme is AnalysisScheme.SHOOTING: + SGM_required_inputs = { + 't_curr': {'units': 's'}, + Dynamic.Mission.DISTANCE: {'units': 'm'}, + } + add_SGM_required_inputs(self, SGM_required_inputs) + self.add_subsystem( 'input_port', VariablesIn(aviary_options=aviary_options), promotes_inputs=['*'], promotes_outputs=['*']) diff --git a/aviary/mission/flops_based/phases/test/test_time_integration_phases.py b/aviary/mission/flops_based/phases/test/test_time_integration_phases.py new file mode 100644 index 000000000..4665096cf --- /dev/null +++ b/aviary/mission/flops_based/phases/test/test_time_integration_phases.py @@ -0,0 +1,204 @@ +import openmdao.api as om +from openmdao.utils.assert_utils import assert_near_equal + +from aviary.interface.methods_for_level2 import AviaryGroup +from aviary.mission.gasp_based.phases.time_integration_traj import FlexibleTraj +from aviary.mission.flops_based.phases.time_integration_phases import \ + SGMHeightEnergy, SGMDetailedTakeoff, SGMDetailedLanding +from aviary.subsystems.premission import CorePreMission +from aviary.utils.functions import set_aviary_initial_values +from aviary.variable_info.enums import Verbosity +from aviary.variable_info.variables import Aircraft, Dynamic, Mission +from aviary.variable_info.variables_in import VariablesIn + +from aviary.interface.default_phase_info.height_energy import aero, prop, geom +from aviary.subsystems.propulsion.engine_deck import EngineDeck +from aviary.utils.process_input_decks import create_vehicle +from aviary.utils.preprocessors import preprocess_propulsion +from aviary.variable_info.variable_meta_data import _MetaData as BaseMetaData + +import warnings +import unittest +import importlib + + +@unittest.skipUnless(importlib.util.find_spec("pyoptsparse") is not None, "pyoptsparse is not installed") +class HE_SGMDescentTestCase(unittest.TestCase): + def setUp(self): + aviary_inputs, initial_guesses = create_vehicle( + 'models/test_aircraft/aircraft_for_bench_FwFm.csv') + aviary_inputs.set_val(Aircraft.Engine.SCALED_SLS_THRUST, val=28690, units="lbf") + aviary_inputs.set_val(Dynamic.Mission.THROTTLE, val=0, units="unitless") + aviary_inputs.set_val(Mission.Takeoff.ROLLING_FRICTION_COEFFICIENT, + val=0.0175, units="unitless") + aviary_inputs.set_val(Mission.Takeoff.BRAKING_FRICTION_COEFFICIENT, + val=0.35, units="unitless") + ode_args = dict(aviary_options=aviary_inputs, core_subsystems=[prop, geom, aero]) + preprocess_propulsion(aviary_inputs, [EngineDeck(options=aviary_inputs)]) + + ode_args['num_nodes'] = 1 + ode_args['subsystem_options'] = {'core_aerodynamics': {'method': 'computed'}} + + self.ode_args = ode_args + self.aviary_inputs = aviary_inputs + self.tol = 1e-5 + + def setup_prob(self, phases) -> om.Problem: + prob = om.Problem() + prob.driver = om.pyOptSparseDriver() + prob.driver.options["optimizer"] = 'IPOPT' + prob.driver.opt_settings['tol'] = 1.0E-6 + prob.driver.opt_settings['mu_init'] = 1e-5 + prob.driver.opt_settings['max_iter'] = 50 + prob.driver.opt_settings['print_level'] = 5 + + aviary_options = self.ode_args['aviary_options'] + subsystems = self.ode_args['core_subsystems'] + + traj = FlexibleTraj( + Phases=phases, + promote_all_auto_ivc=True, + traj_final_state_output=[Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE], + traj_initial_state_input=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + ) + prob.model = AviaryGroup(aviary_options=aviary_options, + aviary_metadata=BaseMetaData) + prob.model.add_subsystem( + 'pre_mission', + CorePreMission(aviary_options=aviary_options, + subsystems=subsystems), + promotes_inputs=['aircraft:*', 'mission:*'], + promotes_outputs=['aircraft:*', 'mission:*'] + ) + prob.model.add_subsystem('traj', traj, + promotes=['aircraft:*', 'mission:*'] + ) + + prob.model.add_subsystem( + "fuel_obj", + om.ExecComp( + "reg_objective = overall_fuel/10000", + reg_objective={"val": 0.0, "units": "unitless"}, + overall_fuel={"units": "lbm"}, + ), + promotes_inputs=[ + ("overall_fuel", Mission.Summary.TOTAL_FUEL_MASS), + ], + promotes_outputs=[("reg_objective", Mission.Objectives.FUEL)], + ) + + prob.model.add_objective(Mission.Objectives.FUEL, ref=1e4) + + prob.model.add_subsystem( + 'input_sink', + VariablesIn(aviary_options=self.aviary_inputs, + meta_data=BaseMetaData), + promotes_inputs=['*'], + promotes_outputs=['*']) + + with warnings.catch_warnings(): + + # Set initial default values for all LEAPS aircraft variables. + set_aviary_initial_values( + prob.model, self.aviary_inputs, meta_data=BaseMetaData) + + warnings.simplefilter("ignore", om.PromotionWarning) + + prob.setup() + + return prob + + def run_simulation(self, phases, initial_values: dict): + prob = self.setup_prob(phases) + + for key, val in initial_values.items(): + prob.set_val(key, **val) + + prob.run_model() + + distance = prob.get_val('traj.distance_final', units='NM')[0] + mass = prob.get_val('traj.mass_final', units='lbm')[0] + alt = prob.get_val('traj.altitude_final', units='ft')[0] + + final_states = {'distance': distance, 'mass': mass, 'altitude': alt} + return final_states + + # def test_takeoff(self): + # initial_values_takeoff = { + # "traj.altitude_initial": {'val': 0, 'units': "ft"}, + # "traj.mass_initial": {'val': 171000, 'units': "lbm"}, + # "traj.distance_initial": {'val': 0, 'units': "NM"}, + # "traj.velocity": {'val': .1, 'units': "m/s"}, + # } + + # ode_args = self.ode_args + # ode_args['friction_key'] = Mission.Takeoff.ROLLING_FRICTION_COEFFICIENT + # brake_release_to_decision = SGMDetailedTakeoff( + # ode_args, + # simupy_args=dict(verbosity=Verbosity.DEBUG,) + # ) + # brake_release_to_decision.clear_triggers() + # brake_release_to_decision.add_trigger(Dynamic.Mission.VELOCITY, value=167.85, units='kn') + + # phases = {'HE': { + # 'ode': brake_release_to_decision, + # 'vals_to_set': {} + # }} + + # final_states = self.run_simulation(phases, initial_values_takeoff) + # # assert_near_equal(final_states['altitude'], 500, self.tol) + # assert_near_equal(final_states['velocity'], 167.85, self.tol) + + def test_cruise(self): + initial_values_cruise = { + "traj.altitude_initial": {'val': 35000, 'units': "ft"}, + "traj.mass_initial": {'val': 171000, 'units': "lbm"}, + "traj.distance_initial": {'val': 0, 'units': "NM"}, + "traj.mach": {'val': .8, 'units': "unitless"}, + } + + SGMCruise = SGMHeightEnergy( + self.ode_args, + phase_name='cruise', + simupy_args=dict(verbosity=Verbosity.QUIET,) + ) + SGMCruise.triggers[0].value = 160000 + + phases = {'HE': { + 'ode': SGMCruise, + 'vals_to_set': {} + }} + + final_states = self.run_simulation(phases, initial_values_cruise) + assert_near_equal(final_states['mass'], 160000, self.tol) + + # def test_landing(self): + # initial_values_landing = { + # "traj.altitude_initial": {'val': 35000, 'units': "ft"}, + # "traj.mass_initial": {'val': 171000, 'units': "lbm"}, + # "traj.distance_initial": {'val': 0, 'units': "NM"}, + # "traj.velocity": {'val': 300, 'units': "m/s"}, + # } + + # ode_args = self.ode_args + # ode_args['friction_key'] = Mission.Takeoff.BRAKING_FRICTION_COEFFICIENT + # phases = {'HE': { + # 'ode': SGMDetailedLanding( + # ode_args, + # simupy_args=dict(verbosity=Verbosity.QUIET,) + # ), + # 'vals_to_set': {} + # }} + + # final_states = self.run_simulation(phases, initial_values_landing) + # assert_near_equal(final_states['altitude'], 0, self.tol) + + +if __name__ == '__main__': + unittest.main() diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py new file mode 100644 index 000000000..7590d469c --- /dev/null +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -0,0 +1,84 @@ +from aviary.mission.flops_based.ode.mission_ODE import MissionODE +from aviary.mission.flops_based.ode.landing_ode import LandingODE +from aviary.mission.flops_based.ode.takeoff_ode import TakeoffODE +from aviary.mission.gasp_based.ode.time_integration_base_classes import SimuPyProblem +from aviary.variable_info.enums import AnalysisScheme +from aviary.variable_info.variables import Dynamic + + +class SGMHeightEnergy(SimuPyProblem): + def __init__( + self, + ode_args, + phase_name='mission', + simupy_args={}, + ): + super().__init__(MissionODE( + analysis_scheme=AnalysisScheme.SHOOTING, + **ode_args), + problem_name=phase_name, + outputs=[], + states=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + alternate_state_rate_names={ + Dynamic.Mission.MASS: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + aviary_options=ode_args['aviary_options'], + **simupy_args) + + self.phase_name = phase_name + self.add_trigger(Dynamic.Mission.MASS, 150000, units='lbm') + + +class SGMDetailedTakeoff(SimuPyProblem): + def __init__( + self, + ode_args, + phase_name='detailed_takeoff', + simupy_args={}, + ): + super().__init__(TakeoffODE( + analysis_scheme=AnalysisScheme.SHOOTING, + **ode_args), + problem_name=phase_name, + outputs=[], + states=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + alternate_state_rate_names={ + Dynamic.Mission.MASS: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + aviary_options=ode_args['aviary_options'], + **simupy_args) + + self.phase_name = phase_name + self.add_trigger(Dynamic.Mission.ALTITUDE, 50, units='ft') + + +class SGMDetailedLanding(SimuPyProblem): + def __init__( + self, + ode_args, + phase_name='detailed_landing', + simupy_args={}, + ): + super().__init__(LandingODE( + analysis_scheme=AnalysisScheme.SHOOTING, + **ode_args), + problem_name=phase_name, + outputs=[], + states=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + alternate_state_rate_names={ + Dynamic.Mission.MASS: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + aviary_options=ode_args['aviary_options'], + **simupy_args) + + self.phase_name = phase_name + self.add_trigger(Dynamic.Mission.ALTITUDE, 0, units='ft') diff --git a/aviary/mission/gasp_based/ode/ascent_eom.py b/aviary/mission/gasp_based/ode/ascent_eom.py index ea6f7c3cd..0cda876c5 100644 --- a/aviary/mission/gasp_based/ode/ascent_eom.py +++ b/aviary/mission/gasp_based/ode/ascent_eom.py @@ -151,6 +151,7 @@ def compute(self, inputs, outputs): * GRAV_ENGLISH_GASP / weight ) + outputs[Dynamic.Mission.FLIGHT_PATH_ANGLE_RATE] = ( (thrust_across_flightpath + incremented_lift - weight * np.cos(gamma)) * GRAV_ENGLISH_GASP diff --git a/aviary/mission/gasp_based/ode/climb_ode.py b/aviary/mission/gasp_based/ode/climb_ode.py index ca9697001..9ca38d361 100644 --- a/aviary/mission/gasp_based/ode/climb_ode.py +++ b/aviary/mission/gasp_based/ode/climb_ode.py @@ -13,8 +13,6 @@ from aviary.subsystems.propulsion.propulsion_builder import PropulsionBuilderBase from aviary.variable_info.enums import AnalysisScheme, AlphaModes, SpeedType from aviary.variable_info.variables import Dynamic -from aviary.mission.ode.specific_energy_rate import SpecificEnergyRate -from aviary.mission.ode.altitude_rate import AltitudeRate class ClimbODE(BaseODE): diff --git a/aviary/mission/gasp_based/ode/descent_ode.py b/aviary/mission/gasp_based/ode/descent_ode.py index fa117f639..ffc498aa3 100644 --- a/aviary/mission/gasp_based/ode/descent_ode.py +++ b/aviary/mission/gasp_based/ode/descent_ode.py @@ -18,8 +18,6 @@ from aviary.variable_info.variables import Dynamic from aviary.subsystems.aerodynamics.aerodynamics_builder import AerodynamicsBuilderBase from aviary.subsystems.propulsion.propulsion_builder import PropulsionBuilderBase -from aviary.mission.ode.specific_energy_rate import SpecificEnergyRate -from aviary.mission.ode.altitude_rate import AltitudeRate class DescentODE(BaseODE): diff --git a/aviary/mission/gasp_based/ode/flight_path_ode.py b/aviary/mission/gasp_based/ode/flight_path_ode.py index 838693eba..3d946c3c7 100644 --- a/aviary/mission/gasp_based/ode/flight_path_ode.py +++ b/aviary/mission/gasp_based/ode/flight_path_ode.py @@ -216,9 +216,6 @@ def setup(self): from aviary.utils.functions import create_printcomp dummy_comp = create_printcomp( all_inputs=[ - Aircraft.Design.OPERATING_MASS, - Aircraft.CrewPayload.PASSENGER_PAYLOAD_MASS, # should be total - Mission.Design.RESERVE_FUEL, Dynamic.Mission.DISTANCE, Dynamic.Mission.THROTTLE, 'required_lift', diff --git a/aviary/mission/gasp_based/ode/time_integration_base_classes.py b/aviary/mission/gasp_based/ode/time_integration_base_classes.py index f8c960629..b0e823843 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -7,23 +7,53 @@ from aviary.mission.gasp_based.ode.params import ParamPort from aviary.variable_info.enums import Verbosity +from aviary.variable_info.variable_meta_data import _MetaData def add_SGM_required_inputs(group: om.Group, inputs_to_add: dict): - blank_execcomp = om.ExecComp() + blank_component = om.ExplicitComponent() for input, details in inputs_to_add.items(): - blank_execcomp.add_input(input, **details) + blank_component.add_input(input, **details) group.add_subsystem('SGM_required_inputs', - blank_execcomp, + blank_component, promotes_inputs=list(inputs_to_add.keys()), ) -# Subproblem used as a basis for forward in time integration phases. +def add_SGM_required_outputs(group: om.Group, outputs_to_add: dict): + iv_comp = om.IndepVarComp() + for output, details in outputs_to_add.items(): + iv_comp.add_output(output, **details) + group.add_subsystem('SGM_required_outputs', + iv_comp, + promotes_outputs=list(outputs_to_add.keys()), + ) + + +class event_trigger(): + def __init__( + self, + state: str, + value: tuple[str, float, int], + units: str, + channel_name: str = None, + ): + self.state = state + self.value = value + self.units = units + if channel_name is None: + channel_name = state + self.channel_name = channel_name + + class SimuPyProblem(SimulationMixin): + # Subproblem used as a basis for forward in time integration phases. def __init__( self, ode, + aviary_options=None, + meta_data=_MetaData, + problem_name='', time_independent=False, t_name="t_curr", states=None, @@ -33,6 +63,7 @@ def __init__( parameters=None, outputs=None, controls=None, + triggers=None, include_state_outputs=False, rate_suffix="_rate", verbosity=Verbosity.QUIET, @@ -64,6 +95,10 @@ def __init__( self.dt = 0.0 prob = om.Problem() + if aviary_options: + from aviary.interface.methods_for_level2 import AviaryGroup + prob.model = AviaryGroup( + aviary_options=aviary_options, aviary_metadata=meta_data) prob.model.add_subsystem( "ODE_group", ode, @@ -75,6 +110,13 @@ def __init__( prob.setup(check=False, force_alloc_complex=True) prob.final_setup() + if triggers is None: + triggers = [] + elif not isinstance(triggers, list): + triggers = [triggers] + self.triggers = triggers + self.event_channel_names = [trigger.channel_name for trigger in triggers] + self.time_independent = time_independent if type(states) is list: states = {state: {'units': None, 'rate': state+rate_suffix, 'rate_units': None} @@ -127,7 +169,7 @@ def __init__( if state_key in states.keys(): # Used to rename an existing state states[val] = states.pop(state_key) else: # Used to add a new state - states[val]: {'units': None, 'rate': val+rate_suffix} + states[val] = {'units': None, 'rate': val+rate_suffix} if alternate_state_rate_names: for state_name, val in alternate_state_rate_names.items(): @@ -160,7 +202,7 @@ def __init__( if outputs is None: outputs = { outp: None - for outp in outputs + for outp in model_outputs if outp not in state_rate_names # + event_names } @@ -172,10 +214,11 @@ def __init__( ).values()))["units"] if include_state_outputs or outputs == {}: # prevent empty outputs - outputs = states + outputs + outputs.update({state: data['units'] for state, data in states.items()}) self.t_name = t_name self.states = states + self.state_names = list(states.keys()) self.parameters = parameters self.outputs = outputs @@ -188,8 +231,11 @@ def __init__( # calls if verbosity.value >= 2: - om.n2(prob, outfile="n2_simupy_problem.html", show_browser=False) - with open('input_list_simupy.txt', 'w') as outfile: + if problem_name: + problem_name = '_'+problem_name + om.n2(prob, outfile="n2_simupy_problem" + + problem_name+".html", show_browser=False) + with open('input_list_simupy'+problem_name+'.txt', 'w') as outfile: prob.model.list_inputs(out_stream=outfile,) print(states) @@ -333,6 +379,43 @@ def update_equation_function(self, t, x, event_channels=None): self.output_nan = True return x + def add_trigger(self, state, value, units=None, channel_name=None): + if units is None: + units = self.states[state]['units'] + elif hasattr(self, units): + units = getattr(self, units) + if channel_name is None: + channel_name = state + + self.triggers.append( + event_trigger(state, value, units, channel_name) + ) + self.event_channel_names.append(channel_name) + self.num_events = len(self.event_channel_names) + + def clear_triggers(self): + self.triggers = [] + self.event_channel_names = [] + self.num_events = 0 + + def event_equation_function(self, t, x): + self.output_equation_function(t, x) + self.compute() + event_values = [self.evaluate_trigger(trigger) for trigger in self.triggers] + # print(event_values) + return np.array(event_values) + + def evaluate_trigger(self, trigger: event_trigger): + trigger_value = trigger.value + if isinstance(trigger_value, str): + if hasattr(self, trigger_value): + trigger_value = getattr(self, trigger_value) + else: + trigger_value = self.get_val( + trigger_value, units=trigger.units).squeeze() + current_value = self.get_val(trigger.state, units=trigger.units).squeeze() + return current_value - trigger_value + @property def get_val(self): return self.prob.get_val @@ -343,12 +426,12 @@ def set_val(self): class SGMTrajBase(om.ExplicitComponent): - def initialize(self): + def initialize(self, verbosity=Verbosity.QUIET): # needs to get passed to each ODE # TODO: param_dict self.options.declare("param_dict", default=ParamPort.param_data) - self.verbosity = Verbosity.QUIET + self.verbosity = verbosity self.max_allowable_time = 1_000_000 self.adjoint_int_opts = DEFAULT_INTEGRATOR_OPTIONS.copy() self.adjoint_int_opts['nsteps'] = 5000 @@ -496,7 +579,7 @@ def compute_traj_loop(self, first_problem, inputs, outputs, t0=0., state0=None): inputs[state_name+"_initial"].squeeze() if state_name in self.traj_initial_state_input else 0. - for state_name in first_problem.states.keys() + for state_name in first_problem.state_names ]).squeeze() while True: @@ -559,7 +642,7 @@ def compute_traj_loop(self, first_problem, inputs, outputs, t0=0., state0=None): outputs[output_name] = sim_results[-1].x[ -1, - list(sim_problems[-1].states.keys()).index(state_name) + sim_problems[-1].state_names.index(state_name) ] for output in self.traj_promote_final_output: @@ -616,13 +699,13 @@ def compute_partials(self, inputs, J): param_deriv = np.zeros(len(param_dict)) if output in self.traj_final_state_output: - costate[list(next_prob.states.keys()).index(output)] = 1. + costate[next_prob.state_names.index(output)] = 1. else: # in self.traj_promote_final_output next_prob.state_equation_function(next_res.t[-1], next_res.x[-1, :]) costate[:] = next_prob.compute_totals( output, - next_prob.states.keys(), + next_prob.state_names, return_format='array' ).squeeze() @@ -663,12 +746,12 @@ def compute_partials(self, inputs, J): num_active_event_channels += 1 dg_dx = np.zeros((1, prob.dim_state)) - if channel_name in prob.states.keys(): - dg_dx[0, list(prob.states.keys()).index(channel_name)] = 1. + if channel_name in prob.state_names: + dg_dx[0, prob.state_names.index(channel_name)] = 1. else: dg_dx[0, :] = prob.compute_totals( [channel_name], - prob.states.keys(), + prob.state_names, return_format='array' ) @@ -704,17 +787,17 @@ def compute_partials(self, inputs, J): # here and co-state assume number of states is only decreasing # forward in time - for state_name in next_prob.states.keys(): - state_idx = list(next_prob.states.keys()).index(state_name) + for state_name in next_prob.state_names: + state_idx = next_prob.state_names.index(state_name) - if state_name in prob.states.keys(): + if state_name in prob.state_names: f_plus[ state_idx - ] = plus_rate[list(prob.states.keys()).index(state_name)] + ] = plus_rate[prob.state_names.index(state_name)] # state_update[ - # list(next_prob.states.keys()).index(state_name) - # ] = x[list(prob.states.keys()).index(state_name)] + # next_prob.state_names.index(state_name) + # ] = x[prob.state_names.index(state_name)] # TODO: make sure index multiplying next_pronb costate # lines up -- since costate is pre-filled to next_prob's @@ -729,7 +812,7 @@ def compute_partials(self, inputs, J): dh_j_dx = prob.compute_totals( [state_name], - prob.states.keys(), + prob.state_names, return_format='array').squeeze() dh_dparam[state_idx, :] = prob.compute_totals( @@ -738,15 +821,15 @@ def compute_partials(self, inputs, J): return_format='array' ).squeeze() - for state_name_2 in prob.states.keys(): + for state_name_2 in prob.state_names: # I'm actually computing dh_dx.T # dh_dx rows are new state, columns are old state # now, dh_dx.T rows are old state, columns are new # so I think this is right dh_dx[ - list(next_prob.states.keys()).index(state_name_2), + next_prob.state_names.index(state_name_2), state_idx, - ] = dh_j_dx[list(prob.states.keys()).index(state_name_2)] + ] = dh_j_dx[prob.state_names.index(state_name_2)] else: state_update[ @@ -760,7 +843,7 @@ def compute_partials(self, inputs, J): state_rate_names = [val['rate'] for _, val in prob.states.items()] df_dx_data[idx, :, :] = prob.compute_totals(state_rate_names, - prob.states.keys(), + prob.state_names, return_format='array').T if param_dict: df_dparam_data[idx, ...] = prob.compute_totals( @@ -875,7 +958,7 @@ def compute_partials(self, inputs, J): # lamda_dot_plus = lamda_dot if self.verbosity is Verbosity.DEBUG: if np.any(state_disc): - print("update is non-zero!", prob, prob.states.keys(), + print("update is non-zero!", prob, prob.state_names, state_disc, costate, lamda_dot) print( "inner product becomes...", @@ -884,7 +967,7 @@ def compute_partials(self, inputs, J): state_disc[None, :] @ dh_dx.T @ lamda_dot_plus[:, None] ) - print("dh_dx for", prob, prob.states.keys(), "\n", dh_dx) + print("dh_dx for", prob, prob.state_names, "\n", dh_dx) print("costate", costate) costate_update_terms = [ dh_dx.T @ costate[:, None], @@ -981,17 +1064,17 @@ def co_state_rate(t, costate, *args): # TODO: do co-states need unit changes? probably not... for state_name in prob.state_names: - costate[list(next_prob.states.keys()).index( - state_name)] = co_res.x[-1, list(prob.states.keys()).index(state_name)] + costate[next_prob.state_names.index( + state_name)] = co_res.x[-1, prob.state_names.index(state_name)] lamda_dot_plus[ - list(next_prob.states.keys()).index(state_name) - ] = lamda_dot_plus_rate[list(prob.states.keys()).index(state_name)] + next_prob.state_names.index(state_name) + ] = lamda_dot_plus_rate[prob.state_names.index(state_name)] for state_to_deriv, metadata in self.traj_initial_state_input.items(): param_name = metadata["name"] J[output_name, param_name] = costate_reses[output][-1].x[ -1, - list(prob.states.keys()).index(state_to_deriv) + prob.state_names.index(state_to_deriv) ] for param_deriv_val, param_deriv_name in zip(param_deriv, param_dict): J[output_name, param_deriv_name] = param_deriv_val diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index 0f6a4962a..7e9841374 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -1,17 +1,14 @@ import numpy as np -import openmdao.api as om from aviary.mission.gasp_based.ode.accel_ode import AccelODE from aviary.mission.gasp_based.ode.ascent_ode import AscentODE from aviary.mission.gasp_based.ode.climb_ode import ClimbODE from aviary.mission.gasp_based.ode.descent_ode import DescentODE from aviary.mission.gasp_based.ode.flight_path_ode import FlightPathODE -from aviary.mission.gasp_based.phases.landing_group import LandingSegment from aviary.mission.gasp_based.ode.groundroll_ode import GroundrollODE from aviary.mission.gasp_based.ode.rotation_ode import RotationODE from aviary.mission.gasp_based.ode.time_integration_base_classes import SimuPyProblem -from aviary.utils.aviary_values import AviaryValues -from aviary.variable_info.enums import AlphaModes, AnalysisScheme, SpeedType, Verbosity +from aviary.variable_info.enums import AlphaModes, AnalysisScheme, SpeedType from aviary.variable_info.variables import Aircraft, Dynamic, Mission @@ -33,6 +30,7 @@ def __init__( super().__init__( GroundrollODE(analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), + problem_name=phase_name, outputs=["normal_force"], states=[ Dynamic.Mission.MASS, @@ -48,14 +46,15 @@ def __init__( self.phase_name = phase_name self.VR_value = VR_value - self.VR_units = VR_units + # self.VR_units = VR_units self.event_channel_names = [Dynamic.Mission.VELOCITY] self.num_events = len(self.event_channel_names) + # self.add_trigger(Dynamic.Mission.VELOCITY, "VR_value", units='ft/s') def event_equation_function(self, t, x): self.time = t self.state = x - return self.get_val("velocity", units='ft/s') - self.VR_value + return self.get_val(Dynamic.Mission.VELOCITY, units='ft/s') - self.VR_value class SGMRotation(SimuPyProblem): @@ -73,11 +72,13 @@ def __init__( ): super().__init__( RotationODE(analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), + problem_name=phase_name, outputs=["normal_force", "alpha"], states=[ Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, + Dynamic.Mission.VELOCITY, ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ @@ -86,15 +87,7 @@ def __init__( ) self.phase_name = phase_name - self.event_channel_names = ["normal_force"] - self.num_events = len(self.event_channel_names) - - def event_equation_function(self, t, x): - self.output_equation_function(t, x) - self.compute() - norm_force = self.get_val("normal_force", units="lbf") + 0.0 - - return norm_force + self.add_trigger("normal_force", 0, units='lbf') # TODO : turn these into parameters? inputs? they need to match between @@ -126,6 +119,7 @@ def __init__( super().__init__( AscentODE(analysis_scheme=AnalysisScheme.SHOOTING, alpha_mode=alpha_mode, **ode_args), + problem_name=phase_name, outputs=[ "load_factor", "fuselage_pitch", @@ -136,6 +130,9 @@ def __init__( Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, + Dynamic.Mission.VELOCITY, + Dynamic.Mission.FLIGHT_PATH_ANGLE, + "alpha", ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ @@ -202,20 +199,21 @@ def __init__( simupy_args={}, ): self.ode_args = ode_args - super().__init__(alpha_mode=AlphaModes.DEFAULT, ode_args=ode_args, simupy_args=simupy_args) + super().__init__(phase_name=phase_name, alpha_mode=AlphaModes.DEFAULT, + ode_args=ode_args, simupy_args=simupy_args) self.phase_name = phase_name self.fuselage_pitch_max = fuselage_pitch_max - ode0 = SGMAscent(alpha_mode=AlphaModes.DEFAULT, + ode0 = SGMAscent(alpha_mode=AlphaModes.DEFAULT, phase_name='ascent_ode0', ode_args=ode_args, simupy_args=simupy_args) - rotation = SGMAscent(alpha_mode=AlphaModes.ROTATION, + rotation = SGMAscent(alpha_mode=AlphaModes.ROTATION, phase_name='ascent_rotation', ode_args=ode_args, simupy_args=simupy_args) - load_factor = SGMAscent(alpha_mode=AlphaModes.LOAD_FACTOR, + load_factor = SGMAscent(alpha_mode=AlphaModes.LOAD_FACTOR, phase_name='ascent_load_factor', ode_args=ode_args, simupy_args=simupy_args) - fuselage_pitch = SGMAscent( - alpha_mode=AlphaModes.FUSELAGE_PITCH, ode_args=ode_args, simupy_args=simupy_args) - decel = SGMAscent(alpha_mode=AlphaModes.DECELERATION, + fuselage_pitch = SGMAscent(alpha_mode=AlphaModes.FUSELAGE_PITCH, phase_name='ascent_pitch', + ode_args=ode_args, simupy_args=simupy_args) + decel = SGMAscent(alpha_mode=AlphaModes.DECELERATION, phase_name='ascent_decel', ode_args=ode_args, simupy_args=simupy_args) self.odes = (ode0, rotation, load_factor, fuselage_pitch, decel,) @@ -374,11 +372,13 @@ def __init__( ode = AccelODE(analysis_scheme=AnalysisScheme.SHOOTING, **ode_args) super().__init__( ode, + problem_name=phase_name, outputs=["EAS", "mach", "alpha"], states=[ Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, + Dynamic.Mission.VELOCITY, ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ @@ -386,16 +386,7 @@ def __init__( **simupy_args, ) self.phase_name = phase_name - self.VC_value = VC_value - self.VC_units = VC_units - self.event_channel_names = [ - "EAS", - ] - self.num_events = len(self.event_channel_names) - - def event_equation_function(self, t, x): - self.output_equation_function(t, x) - return self.get_val("EAS", units=self.VC_units) - self.VC_value + self.add_trigger("EAS", VC_value, units=VC_units) class SGMClimb(SimuPyProblem): @@ -433,6 +424,7 @@ def __init__( self.alt_trigger_units = alt_trigger_units super().__init__( ode, + problem_name=phase_name, outputs=[ "alpha", Dynamic.Mission.FLIGHT_PATH_ANGLE, @@ -456,25 +448,10 @@ def __init__( **simupy_args, ) self.phase_name = phase_name - self.event_channel_names = [ - Dynamic.Mission.ALTITUDE, - self.speed_trigger_name, - ] - self.num_events = len(self.event_channel_names) - - def event_equation_function(self, t, x): - self.output_equation_function(t, x) - alt = self.get_val(Dynamic.Mission.ALTITUDE, - units=self.alt_trigger_units).squeeze() - alt_trigger = self.get_val("alt_trigger", units=self.alt_trigger_units).squeeze() - - speed = self.get_val( - self.speed_trigger_name, units=self.speed_trigger_units - ).squeeze() - speed_trigger = self.get_val( - "speed_trigger", units=self.speed_trigger_units - ).squeeze() - return np.array([alt - alt_trigger, speed - speed_trigger]) + self.add_trigger(Dynamic.Mission.ALTITUDE, "alt_trigger", + units=self.alt_trigger_units) + self.add_trigger(self.speed_trigger_name, "speed_trigger", + units="speed_trigger_units") class SGMCruise(SimuPyProblem): @@ -504,6 +481,7 @@ def __init__( super().__init__( ode, + problem_name=phase_name, outputs=[ "alpha", # ? "lift", @@ -517,6 +495,7 @@ def __init__( Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, + Dynamic.Mission.VELOCITY, ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ @@ -525,26 +504,9 @@ def __init__( ) self.phase_name = phase_name - self.event_channel_names = [ - Dynamic.Mission.MASS, - Dynamic.Mission.DISTANCE, - ] - self.num_events = len(self.event_channel_names) - - def event_equation_function(self, t, x): - self.output_equation_function(t, x) - current_mass = self.get_val(Dynamic.Mission.MASS, units="lbm").squeeze() - mass_trigger = self.get_val('mass_trigger.mass_trigger', units="lbm").squeeze() - - distance = self.get_val(Dynamic.Mission.DISTANCE, - units=self.distance_trigger_units).squeeze() - distance_trigger = self.get_val( - "distance_trigger", units=self.distance_trigger_units).squeeze() - - return np.array([ - current_mass - mass_trigger, - distance - distance_trigger - ]) + self.add_trigger(Dynamic.Mission.MASS, 'mass_trigger.mass_trigger', units='lbm') + self.add_trigger(Dynamic.Mission.DISTANCE, "distance_trigger", + units=self.distance_trigger_units) class SGMDescent(SimuPyProblem): @@ -582,6 +544,7 @@ def __init__( self.alt_trigger_units = alt_trigger_units super().__init__( ode, + problem_name=phase_name, outputs=[ "alpha", "required_lift", @@ -604,22 +567,7 @@ def __init__( ) self.phase_name = phase_name - self.event_channel_names = [ - Dynamic.Mission.ALTITUDE, - self.speed_trigger_name, - ] - self.num_events = len(self.event_channel_names) - - def event_equation_function(self, t, x): - self.output_equation_function(t, x) - alt = self.get_val(Dynamic.Mission.ALTITUDE, - units=self.alt_trigger_units).squeeze() - alt_trigger = self.get_val("alt_trigger", units=self.alt_trigger_units).squeeze() - - speed = self.get_val( - self.speed_trigger_name, units=self.speed_trigger_units - ).squeeze() - speed_trigger = self.get_val( - "speed_trigger", units=self.speed_trigger_units - ).squeeze() - return np.array([alt - alt_trigger, speed - speed_trigger]) + self.add_trigger(Dynamic.Mission.ALTITUDE, "alt_trigger", + units=self.alt_trigger_units) + self.add_trigger(self.speed_trigger_name, "speed_trigger", + units=self.speed_trigger_units) diff --git a/aviary/mission/gasp_based/phases/time_integration_traj.py b/aviary/mission/gasp_based/phases/time_integration_traj.py index a1e4e9dca..25b5c4a5c 100644 --- a/aviary/mission/gasp_based/phases/time_integration_traj.py +++ b/aviary/mission/gasp_based/phases/time_integration_traj.py @@ -1,11 +1,7 @@ import numpy as np import openmdao.api as om -from aviary.mission.gasp_based.ode.time_integration_base_classes import SimuPyProblem, SGMTrajBase -from aviary.mission.gasp_based.phases.time_integration_phases import SGMGroundroll, SGMRotation, SGMAscentCombined, SGMAccel, SGMClimb, SGMCruise, SGMDescent - -from aviary.variable_info.enums import SpeedType -from aviary.utils.aviary_values import AviaryValues -from aviary.variable_info.variables import Aircraft, Mission, Dynamic +from aviary.mission.gasp_based.ode.time_integration_base_classes import SGMTrajBase +from aviary.mission.gasp_based.phases.time_integration_phases import SGMGroundroll, SGMRotation class TimeIntegrationTrajBase(SGMTrajBase): @@ -72,8 +68,9 @@ def compute(self, inputs, outputs): ode_index = 0 sim_gen = self.compute_traj_loop(self.ODEs[0], inputs, outputs) - print('*'*40) - print('Starting: '+self.ODEs[ode_index].phase_name) + if self.verbosity.value >= 1: + print('*'*40) + print('Starting: '+self.ODEs[ode_index].phase_name) for current_problem, sim_result in sim_gen: t_final = sim_result.t[-1] x_final = sim_result.x[-1, :] @@ -86,308 +83,25 @@ def compute(self, inputs, outputs): except IndexError: next_problem = None - print('Finished: '+current_problem.phase_name) + if self.verbosity.value >= 1: + print('Finished: '+current_problem.phase_name) + if next_problem is not None: if type(current_problem) is SGMGroundroll: next_problem.prob.set_val("start_rotation", t_start_rotation) elif type(current_problem) is SGMRotation: next_problem.rotation.set_val("start_rotation", t_start_rotation) - print('Starting: '+next_problem.phase_name) + + if self.verbosity.value >= 1: + print('Starting: '+next_problem.phase_name) sim_gen.send(next_problem) else: - print('Reached the end of the Trajectory!') + if self.verbosity.value >= 1: + print('Reached the end of the Trajectory!') sim_gen.close() break - print('t_final', t_final) - print('x_final', x_final) - print(self.ODEs[-1].states) - -# class SGMTraj1(TimeIntegrationTrajBase): -# ''' -# This combines the phases from brake release to landing -# ''' - -# def initialize(self): -# super().initialize() - -# def setup(self): -# """ -# API requirements: -# pass ODE's, -# next_problem = f(current_problem, current_result) -# initial_state/time/etc -# next_state from last state/output/event information - -# pass in terminal and integrand output functions with derivatives (components) -# -- anything special for final state, final time? -# declare initial state(s) as parameters to take derivative wrt -# assume all other inputs are parameters for deriv? -# """ - -# # distinction is really "meta" inputs to the trajectory itself and ODE -# # parameters that should get automatically passed (previously `passed_inputs`) -# # may also need to just have a specific call out for initial condition -# # parameters? - -# # the "meta" inputs are really initial condition and/or event termination -# # paramters, like initial states or state triggers -# # need an API for that -# # I guess there are similarly at least two categories of outputs, terminal and -# # integral terms. or just say ode must provide it as a state so theyre all meyer -# # terms - -# # actual ODE setup -# groundroll = SGMGroundroll( -# ode_args=(self.options["ode_args_pyc"] if "groundroll" in -# self.options["pyc_phases"] else self.options["ode_args"]) -# ) -# rotation = SGMRotation( -# ode_args=(self.options["ode_args_pyc"] if "groundroll" in -# self.options["pyc_phases"] else self.options["ode_args"]) -# ) -# ascent = SGMAscentCombined( -# ode_args=(self.options["ode_args_pyc"] if "ascent" in -# self.options["pyc_phases"] else self.options["ode_args"]) -# ) -# accel = SGMAccel( -# ode_args=(self.options["ode_args_pyc"] if "accel" in -# self.options["pyc_phases"] else self.options["ode_args"]) -# ) -# climb1 = SGMClimb( -# input_speed_type=SpeedType.EAS, -# input_speed_units="kn", -# ode_args=(self.options["ode_args_pyc"] if "climb" in -# self.options["pyc_phases"] else self.options["ode_args"]) -# ) - -# climb2 = SGMClimb( -# input_speed_type=SpeedType.EAS, -# input_speed_units="kn", -# ode_args=(self.options["ode_args_pyc"] if "climb" in -# self.options["pyc_phases"] else self.options["ode_args"]) -# ) -# climb3 = SGMClimb( -# input_speed_type=SpeedType.MACH, -# input_speed_units="unitless", -# ode_args=(self.options["ode_args_pyc"] if "climb" in -# self.options["pyc_phases"] else self.options["ode_args"]) -# ) -# self.setup_params( -# ODEs=[ -# groundroll, rotation, ascent, accel, climb1, climb2, climb3 -# ], -# traj_final_state_output=[Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE,], -# traj_promote_final_output=[Dynamic.Mission.ALTITUDE_RATE, -# Dynamic.Mission.FLIGHT_PATH_ANGLE, "TAS"], - -# traj_initial_state_input=[Dynamic.Mission.MASS], -# traj_event_trigger_input=[ -# # specify ODE, output_name, with units that SimuPyProblem expects -# # assume event function is of form ODE.output_name - value -# # third key is event_idx associated with input -# (groundroll, "TAS", 0,), -# (climb3, Dynamic.Mission.ALTITUDE, 0,), -# ], -# ) - -# def compute(self, inputs, outputs): -# self.compute_params(inputs) - -# # ODE setup -# (groundroll, rotation, ascent, accel, climb1, climb2, climb3) = self.ODEs - -# cruise_alt = inputs["SGMClimb_"+Dynamic.Mission.ALTITUDE+"_trigger"] -# cruise_mach = self.options["cruise_mach"] - -# groundroll.VR_value = inputs["SGMGroundroll_TAS_trigger"] - -# climb1.set_val("alt_trigger", 10_000, units="ft") -# climb2.set_val("alt_trigger", cruise_alt, units="ft") -# climb3.set_val("alt_trigger", cruise_alt, units="ft") - -# speed_string = {SpeedType.EAS: 'EAS', -# SpeedType.TAS: 'TAS', SpeedType.MACH: 'mach'} -# climb1.set_val(speed_string[climb1.input_speed_type], -# 250., units=climb1.input_speed_units) -# climb2.set_val(speed_string[climb2.input_speed_type], -# 270., units=climb2.input_speed_units) -# climb3.set_val(speed_string[climb3.input_speed_type], cruise_mach, -# units=climb3.input_speed_units) - -# climb1.set_val("speed_trigger", cruise_mach, units=None) -# climb2.set_val("speed_trigger", cruise_mach, units=None) -# climb3.set_val("speed_trigger", 0.0, units=None) - -# for t_var in ["t_init_gear", "t_init_flaps"]: -# ascent.set_val(t_var, 10_000.) -# ascent.rotation.set_val("start_rotation", 10_000.) - -# sim_gen = self.compute_traj_loop(groundroll, inputs, outputs) - -# for current_problem, sim_result in sim_gen: - -# t = sim_result.t[-1] -# x = sim_result.x[-1, :] -# event_idx = np.argmin(np.abs(sim_result.e[-1, :])) - -# # trajectory-specific phase switching -# if current_problem is groundroll: -# if self.verbosity is Verbosity.DEBUG: -# print("starting rotation") -# rotation.prob.set_val("start_rotation", t) -# ascent.rotation.set_val("start_rotation", t) -# next_problem = rotation -# elif current_problem is rotation: -# if self.verbosity is Verbosity.DEBUG: -# print("starting ascent") -# next_problem = ascent -# elif current_problem is ascent: -# next_problem = accel -# if self.verbosity is Verbosity.DEBUG: -# print("starting accel") -# elif current_problem is accel: -# if self.verbosity is Verbosity.DEBUG: -# print("climb1") -# next_problem = climb1 -# elif current_problem is climb1: -# if self.verbosity is Verbosity.DEBUG: -# print("climb2") -# if event_idx != 0: -# if self.verbosity is Verbosity.DEBUG: -# print("expected to hit alt trigger for climb1") -# next_problem = climb2 -# elif current_problem is climb2: -# if self.verbosity is Verbosity.DEBUG: -# print("climb3") -# if event_idx != 1: -# if self.verbosity is Verbosity.DEBUG: -# print("expected to hit speed trigger for climb2") -# next_problem = climb3 -# elif current_problem is climb3: -# if self.verbosity is Verbosity.DEBUG: -# print("climb ending") -# if event_idx != 0: -# if self.verbosity is Verbosity.DEBUG: -# print("expected to hit alt trigger for climb3") -# next_problem = None -# else: -# if self.verbosity is Verbosity.DEBUG: -# print("unexpected termination") -# next_problem = None - -# if next_problem is not None: -# sim_gen.send(next_problem) -# else: -# sim_gen.close() -# break - - -# class SGMTraj2(TimeIntegrationTrajBase): -# ''' -# This combines the phases from end of cruise to landing -# ''' - -# def setup(self): -# """ -# API requirements: -# pass ODE's, -# next_problem = f(current_problem, current_result) -# initial_state/time/etc -# next_state from last state/output/event information - -# pass in terminal and integrand output functions with derivatives (components) -# -- anything special for final state, final time? -# declare initial state(s) as parameters to take derivative wrt -# assume all other inputs are parameters for deriv? -# """ - -# desc1 = SGMDescent( -# input_speed_type=SpeedType.MACH, -# input_speed_units="unitless", -# speed_trigger_units='kn', -# ode_args=self.options["ode_args"], -# ) - -# desc2 = SGMDescent( -# input_speed_type=SpeedType.EAS, -# input_speed_units="kn", -# speed_trigger_units='kn', -# ode_args=self.options["ode_args"], -# ) -# desc3 = SGMDescent( -# input_speed_type=SpeedType.EAS, -# input_speed_units="kn", -# speed_trigger_units='kn', -# ode_args=self.options["ode_args"], -# ) - -# self.setup_params( -# ODEs=[desc1, desc2, desc3], -# traj_final_state_output=[Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE], -# traj_initial_state_input=[Dynamic.Mission.MASS, -# Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE], -# ) -# self.declare_partials(["*"], ["*"],) - -# def compute(self, inputs, outputs): - -# self.compute_params(inputs) - -# # ODE setup -# desc1, desc2, desc3 = self.ODEs - -# cruise_mach = self.options["cruise_mach"] - -# desc1.set_val("alt_trigger", 10_000, units="ft") -# desc2.set_val("alt_trigger", 10_000, units="ft") -# desc3.set_val("alt_trigger", 1_000, units="ft") - -# speed_string = {SpeedType.EAS: 'EAS', -# SpeedType.TAS: 'TAS', SpeedType.MACH: 'mach'} -# desc1.set_val(speed_string[desc1.input_speed_type], cruise_mach, -# units=desc1.input_speed_units) -# desc2.set_val(speed_string[desc2.input_speed_type], -# 350., units=desc2.input_speed_units) -# desc3.set_val(speed_string[desc3.input_speed_type], -# 250., units=desc3.input_speed_units) - -# desc1.set_val("speed_trigger", 350.0, units=desc1.speed_trigger_units) -# desc2.set_val("speed_trigger", 0.0, units=desc2.speed_trigger_units) -# desc3.set_val("speed_trigger", 0.0, units=desc3.speed_trigger_units) - -# # main loop -# sim_gen = self.compute_traj_loop(desc1, inputs, outputs) - -# for current_problem, sim_result in sim_gen: -# t = sim_result.t[-1] -# x = sim_result.x[-1, :] - -# event_idx = np.argmin(np.abs(sim_result.e[-1, :])) - -# # trajectory-specific phase switching -# if current_problem is desc1: -# if self.verbosity is Verbosity.DEBUG: -# print("desc2") -# if event_idx != 1: -# if self.verbosity is Verbosity.DEBUG: -# print("expected to hit speed trigger for desc1") -# sim_gen.send(desc2) -# elif current_problem is desc2: -# if self.verbosity is Verbosity.DEBUG: -# print("desc3") -# if event_idx != 0: -# if self.verbosity is Verbosity.DEBUG: -# print("expected to hit alt trigger for desc2") -# sim_gen.send(desc3) -# elif current_problem is desc3: -# if self.verbosity is Verbosity.DEBUG: -# print("desc ending") -# if event_idx != 0: -# if self.verbosity is Verbosity.DEBUG: -# print("expected to hit alt trigger for desc3") -# sim_gen.close() -# else: -# if self.verbosity is Verbosity.DEBUG: -# print("unexpected termination") -# sim_gen.close() + if self.verbosity.value >= 1: + print('t_final', t_final) + print('x_final', x_final) + print(self.ODEs[-1].states) diff --git a/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py b/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py index db7bde97d..5e6c5191e 100644 --- a/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py +++ b/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py @@ -7,6 +7,7 @@ from aviary.interface.default_phase_info.two_dof import phase_info from aviary.interface.methods_for_level1 import run_aviary from aviary.variable_info.variables import Aircraft, Mission +from aviary.variable_info.enums import AnalysisScheme @use_tempdirs @@ -60,8 +61,30 @@ def test_bench_GwGm_SNOPT(self): assert_near_equal(prob.get_val("traj.desc2.timeseries.distance")[-1], 3675.0, tolerance=rtol) + @require_pyoptsparse(optimizer="IPOPT") + def test_bench_GwGm_shooting(self): + local_phase_info = deepcopy(phase_info) + prob = run_aviary('models/test_aircraft/aircraft_for_bench_GwGm.csv', + local_phase_info, optimizer='IPOPT', run_driver=False, + analysis_scheme=AnalysisScheme.SHOOTING) + + rtol = 0.01 + + assert_near_equal(prob.get_val(Mission.Design.GROSS_MASS), + 174039., tolerance=rtol) + + assert_near_equal(prob.get_val(Aircraft.Design.OPERATING_MASS), + 95509, tolerance=rtol) + + assert_near_equal(prob.get_val('traj.distance_final', units='NM'), + 3774.3, tolerance=rtol) + + assert_near_equal(prob.get_val('traj.mass_final', units='lbm'), + 136823.47, tolerance=rtol) + if __name__ == '__main__': # unittest.main() test = ProblemPhaseTestCase() test.test_bench_GwGm_SNOPT() + test.test_bench_GwGm_shooting()