From d4fd2ac73c55a1f5308fcf52a3cafd84f06e9914 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Thu, 30 Nov 2023 13:57:10 -0800 Subject: [PATCH 01/34] Initial commit for Height Energy SGM phase --- aviary/mission/flops_based/ode/mission_EOM.py | 3 + aviary/mission/flops_based/ode/mission_ODE.py | 47 +++++- .../flops_based/ode/specific_energy_rate.py | 3 + .../phases/time_integration_phases.py | 156 ++++++++++++++++++ .../phases/time_integration_phases.py | 3 +- 5 files changed, 209 insertions(+), 3 deletions(-) create mode 100644 aviary/mission/flops_based/phases/time_integration_phases.py diff --git a/aviary/mission/flops_based/ode/mission_EOM.py b/aviary/mission/flops_based/ode/mission_EOM.py index 946a85468..c679ed093 100644 --- a/aviary/mission/flops_based/ode/mission_EOM.py +++ b/aviary/mission/flops_based/ode/mission_EOM.py @@ -4,12 +4,15 @@ from aviary.mission.flops_based.ode.range_rate import RangeRate from aviary.mission.flops_based.ode.specific_energy_rate import SpecificEnergyRate from aviary.variable_info.variables import Dynamic +from aviary.variable_info.enums import AnalysisScheme class MissionEOM(om.Group): def initialize(self): self.options.declare('num_nodes', types=int, desc='Number of nodes to be evaluated in the RHS') + self.options.declare("analysis_scheme", types=AnalysisScheme, default=AnalysisScheme.COLLOCATION, + desc="The analysis method that will be used to close the trajectory; for example collocation or time integration") def setup(self): nn = self.options['num_nodes'] diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index 0ce3b3541..2a75a2cfd 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -9,8 +9,9 @@ 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): @@ -42,6 +43,12 @@ def initialize(self): self.options.declare( 'use_actual_takeoff_mass', default=False, desc='flag to use actual takeoff mass in the climb phase, otherwise assume 100 kg fuel burn') + 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 @@ -51,6 +58,14 @@ def setup(self): subsystem_options = options['subsystem_options'] engine_count = len(aviary_options.get_val('engine_models')) + blank_execcomp = om.ExecComp() + blank_execcomp.add_input('t_curr') + blank_execcomp.add_input(Dynamic.Mission.RANGE, units='m') + self.add_subsystem('t_curr', + blank_execcomp, + promotes_inputs=['t_curr', Dynamic.Mission.RANGE] + ) + self.add_subsystem( 'input_port', VariablesIn(aviary_options=aviary_options, @@ -124,6 +139,36 @@ def setup(self): Dynamic.Mission.ALTITUDE_RATE, Dynamic.Mission.ALTITUDE_RATE_MAX, Dynamic.Mission.RANGE_RATE]) + self.set_input_defaults(Mission.Design.GROSS_MASS, val=1, units='kg') + self.set_input_defaults( + Aircraft.Fuselage.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.Fuselage.FINENESS, val=1, units='unitless') + self.set_input_defaults(Aircraft.Fuselage.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults( + Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.VerticalTail.FINENESS, val=1, units='unitless') + self.set_input_defaults(Aircraft.VerticalTail.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults( + Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.HorizontalTail.FINENESS, + val=1, units='unitless') + self.set_input_defaults( + Aircraft.HorizontalTail.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults(Aircraft.Wing.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.Wing.FINENESS, val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults( + Aircraft.Wing.SPAN_EFFICIENCY_FACTOR, val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.TAPER_RATIO, val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.THICKNESS_TO_CHORD, + val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.SWEEP, val=1, units='deg') + self.set_input_defaults(Aircraft.Wing.ASPECT_RATIO, val=1, units='unitless') + self.set_input_defaults( + Aircraft.Design.LIFT_DEPENDENT_DRAG_COEFF_FACTOR, val=1, units='unitless') + self.set_input_defaults( + Aircraft.Design.ZERO_LIFT_DRAG_COEFF_FACTOR, val=1, 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') diff --git a/aviary/mission/flops_based/ode/specific_energy_rate.py b/aviary/mission/flops_based/ode/specific_energy_rate.py index 41002d0df..62f0c5cf3 100644 --- a/aviary/mission/flops_based/ode/specific_energy_rate.py +++ b/aviary/mission/flops_based/ode/specific_energy_rate.py @@ -2,6 +2,7 @@ import openmdao.api as om from aviary.constants import GRAV_METRIC_FLOPS as gravity from aviary.variable_info.variables import Dynamic +from aviary.variable_info.enums import AnalysisScheme class SpecificEnergyRate(om.ExplicitComponent): @@ -12,6 +13,8 @@ class SpecificEnergyRate(om.ExplicitComponent): def initialize(self): self.options.declare('num_nodes', types=int) + self.options.declare("analysis_scheme", types=AnalysisScheme, default=AnalysisScheme.COLLOCATION, + desc="The analysis method that will be used to close the trajectory; for example collocation or time integration") def setup(self): nn = self.options['num_nodes'] 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..30c227bc8 --- /dev/null +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -0,0 +1,156 @@ +import numpy as np +import dymos as dm +import openmdao.api as om + +from aviary.mission.flops_based.ode.mission_ODE import MissionODE +from aviary.mission.flops_based.phases.cruise_phase import Cruise +from aviary.mission.gasp_based.phases.time_integration_traj import FlexibleTraj +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 +from aviary.variable_info.variables import Aircraft, Dynamic, Mission + + +class SGMHeightEnergy(SimuPyProblem): + def __init__(self, + ode_args={}, + simupy_args={},): + super().__init__(MissionODE( + analysis_scheme=AnalysisScheme.SHOOTING, + **ode_args), + alternate_state_names={Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS, + 'TAS': 'TAS'}, + alternate_state_rate_names={ + Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + **simupy_args) + + self.event_channel_names = [ + # Dynamic.Mission.DISTANCE, + Dynamic.Mission.MASS, + ] + self.num_events = len(self.event_channel_names) + + def event_equation_function(self, t, x): + self.output_equation_function(t, x) + 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() + + current_mass = self.get_val(Dynamic.Mission.MASS, units="lbm").squeeze() + mass_trigger = 150000 + + return np.array([ + current_mass - mass_trigger + ]) + + +def test_phase(phases): + 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 + + traj = FlexibleTraj( + Phases=phases, + traj_final_state_output=[Dynamic.Mission.MASS, + Dynamic.Mission.RANGE, Dynamic.Mission.ALTITUDE], + traj_initial_state_input=[ + Dynamic.Mission.MASS, + Dynamic.Mission.RANGE, + Dynamic.Mission.ALTITUDE, + ], + ) + prob.model = om.Group() + prob.model.add_subsystem('traj', traj) + + 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.setup() + + prob.run_model() + + +if __name__ == '__main__': + from aviary.interface.default_flops_phases import phase_info + from aviary.subsystems.propulsion.engine_deck import EngineDeck + from aviary.variable_info.variables import Aircraft, Mission, Dynamic + from aviary.utils.UI import create_vehicle + from aviary.utils.preprocessors import preprocess_propulsion + from aviary.variable_info.variable_meta_data import _MetaData as BaseMetaData + + from aviary.subsystems.propulsion.propulsion_builder import CorePropulsionBuilder + from aviary.subsystems.aerodynamics.aerodynamics_builder import CoreAerodynamicsBuilder + core_subsystems = [CorePropulsionBuilder('core_propulsion'), CoreAerodynamicsBuilder( + 'core_aerodynamics', code_origin='FLOPS')] + + aviary_inputs, initial_guesses = create_vehicle( + 'validation_cases/benchmark_tests/bench_4.csv') + aviary_inputs.set_val('debug_mode', False) + aviary_inputs.set_val(Aircraft.Engine.SCALED_SLS_THRUST, val=28690, units="lbf") + aviary_inputs.set_val(Dynamic.Mission.THROTTLE, val=0, units="unitless") + ode_args_tab = dict(aviary_options=aviary_inputs, core_subsystems=core_subsystems) + engine = EngineDeck(options=aviary_inputs) + preprocess_propulsion(aviary_inputs, [engine]) + + # phase_options = phase_info['cruise']['user_options'] + # phase_object = Cruise( + # 'cruise', + # min_altitude=phase_options['min_altitude'], # m + # max_altitude=phase_options['max_altitude'], # m + # min_mach=phase_options['min_mach'], + # max_mach=phase_options['max_mach'], + # fix_initial=phase_options["fix_initial"], + # fix_final=phase_options["fix_final"], + # aero_builder=phase_options["aero_builder"], + # required_altitude_rate=phase_options["required_altitude_rate"], + # mass_f_cruise=phase_options["mass_f_cruise"], + # range_f_cruise=phase_options["range_f_cruise"], + # aviary_options=aviary_inputs, + # ) + # ode_class = MissionODE + # transcription = dm.Radau( + # num_segments=phase_options['num_segments'], order=phase_options['order'], compressed=True) + # external_subsystems = [] + # meta_data = BaseMetaData.copy() + + # phase = phase_object.build_phase( + # ode_class, + # transcription, + # external_subsystems, + # meta_data, + # ) + + ode_args_tab['num_nodes'] = 1 + ode_args_tab['subsystem_options'] = {'core_aerodynamics': {'method': 'computed'}} + phase_kwargs = dict( + ode_args=ode_args_tab, + simupy_args=dict( + DEBUG=False, + blocked_state_names=['engine.nox', 'nox', 'specific_energy'], + ), + ) + phase_vals = { + } + phases = {'HE': { + 'ode': SGMHeightEnergy(**phase_kwargs), + 'vals_to_set': phase_vals + }} + + test_phase(phases) diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index 1d4a784c7..9b619ae42 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -9,8 +9,7 @@ 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 ( - SGMTrajBase, SimuPyProblem) +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 from aviary.variable_info.variables import Aircraft, Dynamic, Mission From ca34560a70b4d8408780e922e6a540dbb2df518b Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 1 Dec 2023 12:17:11 -0800 Subject: [PATCH 02/34] moved empty execcomp into a function to be shared amongst ODEs --- aviary/mission/flops_based/ode/mission_ODE.py | 16 +++---- .../phases/time_integration_phases.py | 46 ++++++------------- .../ode/time_integration_base_classes.py | 19 ++++++-- 3 files changed, 36 insertions(+), 45 deletions(-) diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index 2a75a2cfd..62875b4ef 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -5,7 +5,7 @@ from aviary.mission.flops_based.ode.mission_EOM import MissionEOM from aviary.subsystems.aerodynamics.flops_based.mach_number import MachNumber -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 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 @@ -53,18 +53,18 @@ def initialize(self): 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')) - blank_execcomp = om.ExecComp() - blank_execcomp.add_input('t_curr') - blank_execcomp.add_input(Dynamic.Mission.RANGE, units='m') - self.add_subsystem('t_curr', - blank_execcomp, - promotes_inputs=['t_curr', Dynamic.Mission.RANGE] - ) + if analysis_scheme is AnalysisScheme.SHOOTING: + SGM_required_inputs = { + 't_curr': {'units': 's'}, + Dynamic.Mission.RANGE: {'units': 'm'}, + } + add_SGM_required_inputs(self, SGM_required_inputs) self.add_subsystem( 'input_port', diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index 30c227bc8..b734d1ae1 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -18,8 +18,12 @@ def __init__(self, super().__init__(MissionODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), - alternate_state_names={Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS, - 'TAS': 'TAS'}, + # alternate_state_names={Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS,}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.RANGE, + Dynamic.Mission.ALTITUDE, + ], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args) @@ -57,7 +61,8 @@ def test_phase(phases): traj = FlexibleTraj( Phases=phases, traj_final_state_output=[Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, Dynamic.Mission.ALTITUDE], + Dynamic.Mission.RANGE, + Dynamic.Mission.ALTITUDE], traj_initial_state_input=[ Dynamic.Mission.MASS, Dynamic.Mission.RANGE, @@ -83,6 +88,9 @@ def test_phase(phases): prob.model.add_objective(Mission.Objectives.FUEL, ref=1e4) prob.setup() + prob.set_val("traj.altitude_initial", val=35000, units="ft") + prob.set_val("traj.mass_initial", val=171000, units="lbm") + prob.set_val("traj.range_initial", val=0, units="NM") prob.run_model() @@ -109,41 +117,13 @@ def test_phase(phases): engine = EngineDeck(options=aviary_inputs) preprocess_propulsion(aviary_inputs, [engine]) - # phase_options = phase_info['cruise']['user_options'] - # phase_object = Cruise( - # 'cruise', - # min_altitude=phase_options['min_altitude'], # m - # max_altitude=phase_options['max_altitude'], # m - # min_mach=phase_options['min_mach'], - # max_mach=phase_options['max_mach'], - # fix_initial=phase_options["fix_initial"], - # fix_final=phase_options["fix_final"], - # aero_builder=phase_options["aero_builder"], - # required_altitude_rate=phase_options["required_altitude_rate"], - # mass_f_cruise=phase_options["mass_f_cruise"], - # range_f_cruise=phase_options["range_f_cruise"], - # aviary_options=aviary_inputs, - # ) - # ode_class = MissionODE - # transcription = dm.Radau( - # num_segments=phase_options['num_segments'], order=phase_options['order'], compressed=True) - # external_subsystems = [] - # meta_data = BaseMetaData.copy() - - # phase = phase_object.build_phase( - # ode_class, - # transcription, - # external_subsystems, - # meta_data, - # ) - ode_args_tab['num_nodes'] = 1 ode_args_tab['subsystem_options'] = {'core_aerodynamics': {'method': 'computed'}} phase_kwargs = dict( ode_args=ode_args_tab, simupy_args=dict( - DEBUG=False, - blocked_state_names=['engine.nox', 'nox', 'specific_energy'], + DEBUG=True, + # blocked_state_names=['engine.nox', 'nox', 'specific_energy'], ), ) phase_vals = { 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 d612c1618..3eb5d423e 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -8,6 +8,16 @@ from aviary.mission.gasp_based.ode.params import ParamPort +def add_SGM_required_inputs(group: om.Group, inputs_to_add: dict): + blank_execcomp = om.ExecComp() + for input, details in inputs_to_add.items(): + blank_execcomp.add_input(input, **details) + group.add_subsystem('SGM_required_inputs', + blank_execcomp, + promotes_inputs=list(inputs_to_add.keys()), + ) + + # Subproblem used as a basis for forward in time integration phases. class SimuPyProblem(SimulationMixin): def __init__( @@ -16,11 +26,11 @@ def __init__( time_independent=False, t_name="t_curr", state_names=None, - alternate_state_names=None, + alternate_state_names={}, blocked_state_names=None, state_units=None, state_rate_names=None, - alternate_state_rate_names=None, + alternate_state_rate_names={}, state_rate_units=None, parameter_names=None, parameter_units=None, @@ -97,7 +107,7 @@ def __init__( for name in blocked_state_names: if name in state_names: state_names.remove(name) - if alternate_state_names is not None: + if alternate_state_names: for key, val in alternate_state_names.items(): state_key = key.replace(rate_suffix, "") if state_key in state_names: # Used to rename an existing state @@ -107,7 +117,7 @@ def __init__( if state_rate_names is None: state_rate_names = [state_name + rate_suffix for state_name in state_names] - if alternate_state_names is not None: + if alternate_state_rate_names: state_rate_names = [name if name not in alternate_state_rate_names.keys( ) else alternate_state_rate_names[name] for name in state_rate_names] @@ -215,6 +225,7 @@ def state(self): @state.setter def state(self, value): + print(self.state, value) if np.all(self.state == value): return for state_name, elem_val, unit in zip( From a7c18fc898e31d0db73f5b9d139d9434675c235a Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Mon, 4 Dec 2023 16:45:52 -0800 Subject: [PATCH 03/34] Converted takeoff and landing ODEs to be compatible with SGM --- aviary/mission/flops_based/ode/landing_ode.py | 17 +++++++++++++++++ aviary/mission/flops_based/ode/takeoff_ode.py | 17 +++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/aviary/mission/flops_based/ode/landing_ode.py b/aviary/mission/flops_based/ode/landing_ode.py index c768083bb..555bfb784 100644 --- a/aviary/mission/flops_based/ode/landing_ode.py +++ b/aviary/mission/flops_based/ode/landing_ode.py @@ -12,12 +12,14 @@ 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 from aviary.variable_info.variables import Dynamic as _Dynamic from aviary.variable_info.variables import Mission from aviary.variable_info.variables_in import VariablesIn +from aviary.variable_info.enums import AnalysisScheme Dynamic = _Dynamic.Mission @@ -64,14 +66,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.RANGE: {'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/takeoff_ode.py b/aviary/mission/flops_based/ode/takeoff_ode.py index 73accaa91..b8a378552 100644 --- a/aviary/mission/flops_based/ode/takeoff_ode.py +++ b/aviary/mission/flops_based/ode/takeoff_ode.py @@ -7,12 +7,14 @@ 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 from aviary.variable_info.variables import Dynamic as _Dynamic from aviary.variable_info.variables import Mission from aviary.variable_info.variables_in import VariablesIn +from aviary.variable_info.enums import AnalysisScheme Dynamic = _Dynamic.Mission @@ -62,14 +64,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.RANGE: {'units': 'm'}, + } + add_SGM_required_inputs(self, SGM_required_inputs) + self.add_subsystem( 'input_port', VariablesIn(aviary_options=aviary_options), promotes_inputs=['*'], promotes_outputs=['*']) From 27d6c128ab0acd1159ec29aac69ca4538a17dced Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Wed, 6 Dec 2023 09:08:26 -0800 Subject: [PATCH 04/34] Updated variable pass through from outer problem to SimuPy problem Added required premission susbsytems for testing --- .../phases/time_integration_phases.py | 66 +++++++++++++++++-- .../ode/time_integration_base_classes.py | 19 ++++-- .../phases/time_integration_traj.py | 2 + 3 files changed, 74 insertions(+), 13 deletions(-) diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index b734d1ae1..f80290220 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -9,6 +9,10 @@ from aviary.utils.aviary_values import AviaryValues from aviary.variable_info.enums import AlphaModes, AnalysisScheme, SpeedType from aviary.variable_info.variables import Aircraft, Dynamic, Mission +from aviary.variable_info.variables_in import VariablesIn +from aviary.subsystems.premission import CorePreMission +from aviary.utils.functions import set_aviary_initial_values +import warnings class SGMHeightEnergy(SimuPyProblem): @@ -49,7 +53,7 @@ def event_equation_function(self, t, x): ]) -def test_phase(phases): +def test_phase(phases, ode_args_tab): prob = om.Problem() prob.driver = om.pyOptSparseDriver() prob.driver.options["optimizer"] = 'IPOPT' @@ -58,8 +62,18 @@ def test_phase(phases): prob.driver.opt_settings['max_iter'] = 50 prob.driver.opt_settings['print_level'] = 5 + aviary_options = ode_args_tab['aviary_options'] + subsystems = ode_args_tab['core_subsystems'] + traj = FlexibleTraj( Phases=phases, + traj_promote_initial_input={ + Aircraft.Wing.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + Aircraft.VerticalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + Aircraft.Fuselage.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + Aircraft.Nacelle.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + }, traj_final_state_output=[Dynamic.Mission.MASS, Dynamic.Mission.RANGE, Dynamic.Mission.ALTITUDE], @@ -70,7 +84,20 @@ def test_phase(phases): ], ) prob.model = om.Group() - prob.model.add_subsystem('traj', traj) + 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.promotes('traj', inputs=[ + Aircraft.Wing.CHARACTERISTIC_LENGTH, + Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, + Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, + Aircraft.Fuselage.CHARACTERISTIC_LENGTH, + Aircraft.Nacelle.CHARACTERISTIC_LENGTH,]) prob.model.add_subsystem( "fuel_obj", @@ -87,12 +114,39 @@ def test_phase(phases): prob.model.add_objective(Mission.Objectives.FUEL, ref=1e4) + prob.model.add_subsystem( + 'input_sink', + VariablesIn(aviary_options=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, aviary_inputs, meta_data=BaseMetaData) + + warnings.simplefilter("ignore", om.PromotionWarning) + prob.setup() prob.set_val("traj.altitude_initial", val=35000, units="ft") prob.set_val("traj.mass_initial", val=171000, units="lbm") prob.set_val("traj.range_initial", val=0, units="NM") + # try: prob.run_model() + # except: + # prob.final_setup() + + with open('input_list.txt', 'w') as outfile: + prob.model.list_inputs(out_stream=outfile,) + with open('output_list.txt', 'w') as outfile: + prob.model.list_outputs(out_stream=outfile) + + final_range = prob.get_val('traj.range_final', units='NM')[0] + final_mass = prob.get_val('traj.mass_final', units='lbm')[0] + print(final_range, final_mass) if __name__ == '__main__': @@ -103,10 +157,8 @@ def test_phase(phases): from aviary.utils.preprocessors import preprocess_propulsion from aviary.variable_info.variable_meta_data import _MetaData as BaseMetaData - from aviary.subsystems.propulsion.propulsion_builder import CorePropulsionBuilder - from aviary.subsystems.aerodynamics.aerodynamics_builder import CoreAerodynamicsBuilder - core_subsystems = [CorePropulsionBuilder('core_propulsion'), CoreAerodynamicsBuilder( - 'core_aerodynamics', code_origin='FLOPS')] + from aviary.interface.default_flops_phases import aero, prop, geom + core_subsystems = [prop, geom, aero] aviary_inputs, initial_guesses = create_vehicle( 'validation_cases/benchmark_tests/bench_4.csv') @@ -133,4 +185,4 @@ def test_phase(phases): 'vals_to_set': phase_vals }} - test_phase(phases) + test_phase(phases, ode_args_tab) 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 3eb5d423e..4b6f7b598 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -199,6 +199,8 @@ def __init__( if DEBUG: om.n2(prob, outfile="n2_simupy_problem.html", show_browser=False) + with open('input_list_simupy.txt', 'w') as outfile: + prob.model.list_inputs(out_stream=outfile,) print(state_names) print(self.state_units) print(state_rate_names) @@ -373,6 +375,7 @@ def setup_params( ODEs, traj_final_state_output=None, traj_promote_final_output=None, + traj_promote_initial_input=None, traj_initial_state_input=None, traj_event_trigger_input=None, ): @@ -392,13 +395,18 @@ def setup_params( traj_final_state_output = [] if traj_promote_final_output is None: traj_promote_final_output = [] + if traj_promote_initial_input is None: + traj_promote_initial_input = {} if traj_initial_state_input is None: traj_initial_state_input = [] if traj_event_trigger_input is None: traj_event_trigger_input = [] - for name, kwargs in self.options["param_dict"].items(): + self.traj_promote_initial_input = { + **self.options["param_dict"], **traj_promote_initial_input} + for name, kwargs in self.traj_promote_initial_input.items(): self.add_input(name, **kwargs) + final_suffix = "_final" self.traj_final_state_output = { final_state_output: { @@ -474,17 +482,16 @@ def setup_params( self.declare_partials(["*"], ["*"],) def compute_params(self, inputs): - # parameter pass-through setup - for param_input in self.options["param_dict"].keys(): + for input in self.traj_promote_initial_input.keys(): for ode in self.ODEs: try: - ode.set_val(param_input, inputs[param_input]) + ode.set_val(input, inputs[input]) except KeyError: if self.DEBUG: print( - "*** ParamPort input not found:", + "*** Input not found:", ode, - param_input + input ) pass diff --git a/aviary/mission/gasp_based/phases/time_integration_traj.py b/aviary/mission/gasp_based/phases/time_integration_traj.py index cc136d943..12f0b6791 100644 --- a/aviary/mission/gasp_based/phases/time_integration_traj.py +++ b/aviary/mission/gasp_based/phases/time_integration_traj.py @@ -31,6 +31,7 @@ def initialize(self): self.options.declare('Phases', default=None) self.options.declare('traj_final_state_output', default=None) self.options.declare('traj_promote_final_output', default=None) + self.options.declare('traj_promote_initial_input', default=None) self.options.declare('traj_initial_state_input', default=None) self.options.declare('traj_event_trigger_input', default=None) @@ -47,6 +48,7 @@ def setup(self): traj_final_state_output=self.options['traj_final_state_output'], traj_promote_final_output=self.options['traj_promote_final_output'], + traj_promote_initial_input=self.options['traj_promote_initial_input'], traj_initial_state_input=self.options['traj_initial_state_input'], traj_event_trigger_input=self.options['traj_event_trigger_input'], ) From f3b8fcdad3cead295849e8b0e907ad46b5360ed6 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Wed, 6 Dec 2023 09:58:10 -0800 Subject: [PATCH 05/34] temp for debugging --- aviary/mission/flops_based/ode/mission_ODE.py | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index 62875b4ef..b48ad3dd5 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -197,3 +197,29 @@ def setup(self): ('mass', Dynamic.Mission.MASS) ], promotes_outputs=['initial_mass_residual']) + + if analysis_scheme is AnalysisScheme.SHOOTING: + 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.MASS, + Dynamic.Mission.RANGE, + Dynamic.Mission.THROTTLE, + 'required_lift', + 'load_factor', + Dynamic.Mission.ALTITUDE, + Dynamic.Mission.FLIGHT_PATH_ANGLE, + ], + input_units={ + 'required_lift': 'lbf', + Dynamic.Mission.FLIGHT_PATH_ANGLE: 'deg', + }) + self.add_subsystem( + "dummy_comp", + dummy_comp(), + promotes_inputs=["*"],) + self.set_input_defaults( + Dynamic.Mission.RANGE, val=0, units='NM') From 7d3b46ef7da573c92942f6de3e6c53b5dff51898 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Thu, 7 Dec 2023 16:28:52 -0800 Subject: [PATCH 06/34] temp commit for debugging --- aviary/mission/flops_based/ode/range_rate.py | 1 + 1 file changed, 1 insertion(+) diff --git a/aviary/mission/flops_based/ode/range_rate.py b/aviary/mission/flops_based/ode/range_rate.py index 9ef5d8fca..7cf4b5693 100644 --- a/aviary/mission/flops_based/ode/range_rate.py +++ b/aviary/mission/flops_based/ode/range_rate.py @@ -32,6 +32,7 @@ def compute(self, inputs, outputs): velocity = inputs[Dynamic.Mission.VELOCITY] climb_rate_2 = climb_rate**2 velocity_2 = velocity**2 + print(climb_rate, velocity) if (climb_rate_2 >= velocity_2).any(): raise om.AnalysisError( "WARNING: climb rate exceeds velocity (range_rate.py)") From 2c4701f19780cc289eb198b47291462fc4eb4d5e Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Thu, 28 Dec 2023 17:34:19 -0800 Subject: [PATCH 07/34] Reduced usage of blocked_state_names by specifying states in more phases --- aviary/examples/level2_shooting_traj.py | 1 - .../interface/default_phase_info/gasp_fiti.py | 10 ---- .../phases/time_integration_phases.py | 7 +-- .../phases/time_integration_phases.py | 57 ++++++++++++++----- 4 files changed, 45 insertions(+), 30 deletions(-) diff --git a/aviary/examples/level2_shooting_traj.py b/aviary/examples/level2_shooting_traj.py index 960db5926..abbf6548a 100644 --- a/aviary/examples/level2_shooting_traj.py +++ b/aviary/examples/level2_shooting_traj.py @@ -74,7 +74,6 @@ def run_aviary(aircraft_filename, phase_info, mission_method, mass_method, optim ode_args=prob.ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine.nox', 'nox'], ), ) descent1_vals = { diff --git a/aviary/interface/default_phase_info/gasp_fiti.py b/aviary/interface/default_phase_info/gasp_fiti.py index 88d3a16d6..b2aff1a89 100644 --- a/aviary/interface/default_phase_info/gasp_fiti.py +++ b/aviary/interface/default_phase_info/gasp_fiti.py @@ -17,7 +17,6 @@ def create_gasp_based_ascent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) groundroll_vals = { @@ -29,7 +28,6 @@ def create_gasp_based_ascent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) rotation_vals = {} @@ -38,7 +36,6 @@ def create_gasp_based_ascent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) ascent_vals = { @@ -53,7 +50,6 @@ def create_gasp_based_ascent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) accel_vals = {} @@ -64,7 +60,6 @@ def create_gasp_based_ascent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) climb1_vals = { @@ -79,7 +74,6 @@ def create_gasp_based_ascent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) climb2_vals = { @@ -94,7 +88,6 @@ def create_gasp_based_ascent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) climb3_vals = { @@ -150,7 +143,6 @@ def create_gasp_based_descent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) descent1_vals = { @@ -166,7 +158,6 @@ def create_gasp_based_descent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) descent2_vals = { @@ -182,7 +173,6 @@ def create_gasp_based_descent_phases( ode_args=ode_args, simupy_args=dict( DEBUG=False, - blocked_state_names=['engine_deck.nox', 'nox'], ), ) descent3_vals = { diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index f80290220..3c9dea9cd 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -22,7 +22,6 @@ def __init__(self, super().__init__(MissionODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), - # alternate_state_names={Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS,}, state_names=[ Dynamic.Mission.MASS, Dynamic.Mission.RANGE, @@ -150,14 +149,13 @@ def test_phase(phases, ode_args_tab): if __name__ == '__main__': - from aviary.interface.default_flops_phases import phase_info + from aviary.interface.default_phase_info.flops import phase_info, aero, prop, geom from aviary.subsystems.propulsion.engine_deck import EngineDeck from aviary.variable_info.variables import Aircraft, Mission, Dynamic - from aviary.utils.UI import create_vehicle + 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 - from aviary.interface.default_flops_phases import aero, prop, geom core_subsystems = [prop, geom, aero] aviary_inputs, initial_guesses = create_vehicle( @@ -175,7 +173,6 @@ def test_phase(phases, ode_args_tab): ode_args=ode_args_tab, simupy_args=dict( DEBUG=True, - # blocked_state_names=['engine.nox', 'nox', 'specific_energy'], ), ) phase_vals = { diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index 9b619ae42..7b5865192 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -36,8 +36,13 @@ def __init__( super().__init__( GroundrollODE(analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), output_names=["normal_force"], - alternate_state_names={Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS, - 'TAS': 'TAS'}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + 'TAS' + ], + # state_units=['lbm','nmi','ft','ft/s'], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args, @@ -72,8 +77,12 @@ def __init__( super().__init__( RotationODE(analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), output_names=["normal_force", "alpha"], - alternate_state_names={ - Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args, @@ -126,8 +135,12 @@ def __init__( "normal_force", "alpha", ], - alternate_state_names={ - Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, control_names=control_names, @@ -365,8 +378,12 @@ def __init__( super().__init__( ode, output_names=["EAS", "mach", "alpha"], - alternate_state_names={ - Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args, @@ -431,8 +448,12 @@ def __init__( "drag", Dynamic.Mission.ALTITUDE_RATE, ], - alternate_state_names={ - Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args, @@ -495,8 +516,12 @@ def __init__( "drag", Dynamic.Mission.ALTITUDE_RATE, ], - alternate_state_names={ - Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args, @@ -570,8 +595,12 @@ def __init__( "drag", Dynamic.Mission.ALTITUDE_RATE, ], - alternate_state_names={ - Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL: Dynamic.Mission.MASS}, + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.DISTANCE, + Dynamic.Mission.ALTITUDE, + ], + # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args, From 22856552acb51112ea31760a121452fe5ac8151d Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 29 Dec 2023 13:12:24 -0800 Subject: [PATCH 08/34] Added ability to auto-promote all auto_ivs variables --- .../phases/time_integration_phases.py | 21 +++++++++++-------- .../ode/time_integration_base_classes.py | 16 ++++++++++++++ .../phases/time_integration_traj.py | 2 ++ 3 files changed, 30 insertions(+), 9 deletions(-) diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index 3c9dea9cd..8cd771e22 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -66,13 +66,14 @@ def test_phase(phases, ode_args_tab): traj = FlexibleTraj( Phases=phases, - traj_promote_initial_input={ - Aircraft.Wing.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - Aircraft.VerticalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - Aircraft.Fuselage.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - Aircraft.Nacelle.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - }, + promote_all_auto_ivc=True, + # traj_promote_initial_input={ + # Aircraft.Wing.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + # Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + # Aircraft.VerticalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + # Aircraft.Fuselage.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + # Aircraft.Nacelle.CHARACTERISTIC_LENGTH: {'units': 'ft'}, + # }, traj_final_state_output=[Dynamic.Mission.MASS, Dynamic.Mission.RANGE, Dynamic.Mission.ALTITUDE], @@ -128,10 +129,12 @@ def test_phase(phases, ode_args_tab): warnings.simplefilter("ignore", om.PromotionWarning) - prob.setup() + prob.setup() prob.set_val("traj.altitude_initial", val=35000, units="ft") prob.set_val("traj.mass_initial", val=171000, units="lbm") prob.set_val("traj.range_initial", val=0, units="NM") + prob.set_val("traj.velocity", val=472, units="kn") + prob.set_val("traj.velocity_rate", val=0, units="m/s**2") # try: prob.run_model() @@ -159,7 +162,7 @@ def test_phase(phases, ode_args_tab): core_subsystems = [prop, geom, aero] aviary_inputs, initial_guesses = create_vehicle( - 'validation_cases/benchmark_tests/bench_4.csv') + 'models/test_aircraft/aircraft_for_bench_FwFm.csv') aviary_inputs.set_val('debug_mode', False) aviary_inputs.set_val(Aircraft.Engine.SCALED_SLS_THRUST, val=28690, units="lbf") aviary_inputs.set_val(Dynamic.Mission.THROTTLE, val=0, units="unitless") 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 4b6f7b598..0a923ae40 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -227,6 +227,7 @@ def state(self): @state.setter def state(self, value): + print('state.setter: self.state, value') print(self.state, value) if np.all(self.state == value): return @@ -373,6 +374,7 @@ def initialize(self): def setup_params( self, ODEs, + promote_all_auto_ivc=False, traj_final_state_output=None, traj_promote_final_output=None, traj_promote_initial_input=None, @@ -402,6 +404,20 @@ def setup_params( if traj_event_trigger_input is None: traj_event_trigger_input = [] + if promote_all_auto_ivc: + for ode in ODEs: + ode: SimuPyProblem + ode.auto_ivc_vars = ode.prob.model.list_inputs( + is_indep_var=True, units=True, out_stream=None) + for abs_name, data in ode.auto_ivc_vars: + prom_name = data.pop('prom_name') + if '.' in prom_name: + continue + traj_promote_initial_input[prom_name] = data + # print(traj_promote_initial_input) + # print(len(traj_promote_initial_input)) + # exit() + self.traj_promote_initial_input = { **self.options["param_dict"], **traj_promote_initial_input} for name, kwargs in self.traj_promote_initial_input.items(): diff --git a/aviary/mission/gasp_based/phases/time_integration_traj.py b/aviary/mission/gasp_based/phases/time_integration_traj.py index 12f0b6791..d566635d3 100644 --- a/aviary/mission/gasp_based/phases/time_integration_traj.py +++ b/aviary/mission/gasp_based/phases/time_integration_traj.py @@ -29,6 +29,7 @@ class FlexibleTraj(TimeIntegrationTrajBase): def initialize(self): super().initialize() self.options.declare('Phases', default=None) + self.options.declare('promote_all_auto_ivc', default=False) self.options.declare('traj_final_state_output', default=None) self.options.declare('traj_promote_final_output', default=None) self.options.declare('traj_promote_initial_input', default=None) @@ -45,6 +46,7 @@ def setup(self): self.setup_params( ODEs=ODEs, + promote_all_auto_ivc=self.options['promote_all_auto_ivc'], traj_final_state_output=self.options['traj_final_state_output'], traj_promote_final_output=self.options['traj_promote_final_output'], From 5fc9c65905cc407c28ffcd01ed1c5fc3a13bcd9b Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 29 Dec 2023 13:14:24 -0800 Subject: [PATCH 09/34] Added debug prints --- aviary/mission/flops_based/ode/altitude_rate.py | 3 +++ aviary/mission/flops_based/ode/range_rate.py | 2 ++ aviary/mission/flops_based/ode/specific_energy_rate.py | 3 +++ aviary/subsystems/aerodynamics/flops_based/drag.py | 3 +++ 4 files changed, 11 insertions(+) diff --git a/aviary/mission/flops_based/ode/altitude_rate.py b/aviary/mission/flops_based/ode/altitude_rate.py index 30c045c1d..fa52cd7d2 100644 --- a/aviary/mission/flops_based/ode/altitude_rate.py +++ b/aviary/mission/flops_based/ode/altitude_rate.py @@ -34,6 +34,9 @@ def compute(self, inputs, outputs): acceleration = inputs[Dynamic.Mission.VELOCITY_RATE] velocity = inputs[Dynamic.Mission.VELOCITY] + print('alt rate:') + print('specific_power, acceleration, velocity') + print(specific_power, acceleration, velocity) outputs[Dynamic.Mission.ALTITUDE_RATE] = \ specific_power - (velocity * acceleration) / gravity diff --git a/aviary/mission/flops_based/ode/range_rate.py b/aviary/mission/flops_based/ode/range_rate.py index 7cf4b5693..f69e419c2 100644 --- a/aviary/mission/flops_based/ode/range_rate.py +++ b/aviary/mission/flops_based/ode/range_rate.py @@ -32,6 +32,8 @@ def compute(self, inputs, outputs): velocity = inputs[Dynamic.Mission.VELOCITY] climb_rate_2 = climb_rate**2 velocity_2 = velocity**2 + print('range rate:') + print('climb_rate, velocity') print(climb_rate, velocity) if (climb_rate_2 >= velocity_2).any(): raise om.AnalysisError( diff --git a/aviary/mission/flops_based/ode/specific_energy_rate.py b/aviary/mission/flops_based/ode/specific_energy_rate.py index 62f0c5cf3..658bb853d 100644 --- a/aviary/mission/flops_based/ode/specific_energy_rate.py +++ b/aviary/mission/flops_based/ode/specific_energy_rate.py @@ -46,6 +46,9 @@ def compute(self, inputs, outputs): weight = inputs[Dynamic.Mission.MASS] * gravity outputs[Dynamic.Mission.SPECIFIC_ENERGY_RATE] = velocity * \ (thrust - drag) / weight + print('spec energy rate:') + print('velocity, thrust, drag, weight') + print(velocity, thrust, drag, weight) def setup_partials(self): arange = np.arange(self.options['num_nodes']) diff --git a/aviary/subsystems/aerodynamics/flops_based/drag.py b/aviary/subsystems/aerodynamics/flops_based/drag.py index 0a1c4ee73..75a1b06c1 100644 --- a/aviary/subsystems/aerodynamics/flops_based/drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/drag.py @@ -66,6 +66,9 @@ def compute(self, inputs, outputs): q = inputs[Dynamic.Mission.DYNAMIC_PRESSURE] CD = inputs['drag_coefficient'] + print('simple drag') + print('S,FCDSUB,FCDSUP') + print(S, FCDSUB, FCDSUP) idx_sup = np.where(M >= 1.0) CD_scaled = CD * FCDSUB CD_scaled[idx_sup] = CD[idx_sup] * FCDSUP From 4676eab4e41198176bff77f81cfb93f42572c39b Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 29 Dec 2023 13:35:40 -0800 Subject: [PATCH 10/34] added more prints --- .../aerodynamics/flops_based/compressibility_drag.py | 3 +++ aviary/subsystems/aerodynamics/flops_based/drag.py | 4 ++-- aviary/subsystems/aerodynamics/flops_based/induced_drag.py | 3 +++ .../aerodynamics/flops_based/lift_dependent_drag.py | 3 +++ .../subsystems/aerodynamics/flops_based/skin_friction_drag.py | 3 +++ 5 files changed, 14 insertions(+), 2 deletions(-) diff --git a/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py b/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py index 87d2c6379..825efac2c 100644 --- a/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py @@ -76,6 +76,9 @@ def compute(self, inputs, outputs): outputs['compress_drag_coeff'][idx_super] = cdc_super outputs['compress_drag_coeff'][idx_sub] = cdc_sub + print('compressibility drag:') + print('cdc_sub, cdc_super') + print(cdc_sub, cdc_super) def _compute_supersonic(self, inputs, outputs, idx): """ diff --git a/aviary/subsystems/aerodynamics/flops_based/drag.py b/aviary/subsystems/aerodynamics/flops_based/drag.py index 75a1b06c1..8cea0c931 100644 --- a/aviary/subsystems/aerodynamics/flops_based/drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/drag.py @@ -67,8 +67,8 @@ def compute(self, inputs, outputs): CD = inputs['drag_coefficient'] print('simple drag') - print('S,FCDSUB,FCDSUP') - print(S, FCDSUB, FCDSUP) + print('M, q, CD') + print(M, q, CD) idx_sup = np.where(M >= 1.0) CD_scaled = CD * FCDSUB CD_scaled[idx_sup] = CD[idx_sup] * FCDSUP diff --git a/aviary/subsystems/aerodynamics/flops_based/induced_drag.py b/aviary/subsystems/aerodynamics/flops_based/induced_drag.py index 8ac3f4726..e61278021 100644 --- a/aviary/subsystems/aerodynamics/flops_based/induced_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/induced_drag.py @@ -73,6 +73,9 @@ def compute(self, inputs, outputs): CL = 2.0 * lift / (Sref * gamma * P * mach ** 2) + print('induced drag:') + print('gamma, mach, lift, P, Sref, AR, span_efficiency_factor, SW25, TR') + print(gamma, mach, lift, P, Sref, AR, span_efficiency_factor, SW25, TR) redux, _ = aviary_options.get_item( Aircraft.Wing.SPAN_EFFICIENCY_REDUCTION, (False, None)) diff --git a/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py b/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py index adba49d48..ab251440e 100644 --- a/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py @@ -115,6 +115,9 @@ def compute(self, inputs, outputs): gamma = self.options['gamma'] mach, lift, P, CLDES, MDES, Sref, AR, CAM, SW25, TC = inputs.values() + print('lift dependent drag:') + print('gamma, mach, lift, P, CLDES, MDES, Sref, AR, CAM, SW25, TC') + print(gamma, mach, lift, P, CLDES, MDES, Sref, AR, CAM, SW25, TC) FCDP = np.empty(nn, dtype=mach.dtype) DCDP = np.empty(nn, dtype=mach.dtype) dFCDP_dDELM = np.empty(nn, dtype=mach.dtype) diff --git a/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py b/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py index c7a9edd0f..fb2ff3639 100644 --- a/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py @@ -108,6 +108,9 @@ def compute(self, inputs, outputs): mission_wing_area = inputs[Aircraft.Wing.AREA] laminar_flow = np.any(lam_up > 0.0) or np.any(lam_low > 0.0) + print('skin friction:') + print('cf, Re, fineness, wetted_area, lam_up, lam_low, laminar_flow') + print(cf, Re, fineness, wetted_area, lam_up, lam_low, laminar_flow) if laminar_flow: laminar_upper = _calc_laminar_flow(lam_up) From 6a6f5a611e7c63773d307be83117f5e2a415ecdd Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 29 Dec 2023 14:16:44 -0800 Subject: [PATCH 11/34] Remove debug prints This reverts commit 5fc9c65905cc407c28ffcd01ed1c5fc3a13bcd9b. --- aviary/mission/flops_based/ode/altitude_rate.py | 3 --- aviary/mission/flops_based/ode/range_rate.py | 2 -- aviary/mission/flops_based/ode/specific_energy_rate.py | 3 --- .../aerodynamics/flops_based/compressibility_drag.py | 3 --- aviary/subsystems/aerodynamics/flops_based/drag.py | 3 --- aviary/subsystems/aerodynamics/flops_based/induced_drag.py | 3 --- .../subsystems/aerodynamics/flops_based/lift_dependent_drag.py | 3 --- .../subsystems/aerodynamics/flops_based/skin_friction_drag.py | 3 --- 8 files changed, 23 deletions(-) diff --git a/aviary/mission/flops_based/ode/altitude_rate.py b/aviary/mission/flops_based/ode/altitude_rate.py index fa52cd7d2..30c045c1d 100644 --- a/aviary/mission/flops_based/ode/altitude_rate.py +++ b/aviary/mission/flops_based/ode/altitude_rate.py @@ -34,9 +34,6 @@ def compute(self, inputs, outputs): acceleration = inputs[Dynamic.Mission.VELOCITY_RATE] velocity = inputs[Dynamic.Mission.VELOCITY] - print('alt rate:') - print('specific_power, acceleration, velocity') - print(specific_power, acceleration, velocity) outputs[Dynamic.Mission.ALTITUDE_RATE] = \ specific_power - (velocity * acceleration) / gravity diff --git a/aviary/mission/flops_based/ode/range_rate.py b/aviary/mission/flops_based/ode/range_rate.py index f69e419c2..7cf4b5693 100644 --- a/aviary/mission/flops_based/ode/range_rate.py +++ b/aviary/mission/flops_based/ode/range_rate.py @@ -32,8 +32,6 @@ def compute(self, inputs, outputs): velocity = inputs[Dynamic.Mission.VELOCITY] climb_rate_2 = climb_rate**2 velocity_2 = velocity**2 - print('range rate:') - print('climb_rate, velocity') print(climb_rate, velocity) if (climb_rate_2 >= velocity_2).any(): raise om.AnalysisError( diff --git a/aviary/mission/flops_based/ode/specific_energy_rate.py b/aviary/mission/flops_based/ode/specific_energy_rate.py index 658bb853d..62f0c5cf3 100644 --- a/aviary/mission/flops_based/ode/specific_energy_rate.py +++ b/aviary/mission/flops_based/ode/specific_energy_rate.py @@ -46,9 +46,6 @@ def compute(self, inputs, outputs): weight = inputs[Dynamic.Mission.MASS] * gravity outputs[Dynamic.Mission.SPECIFIC_ENERGY_RATE] = velocity * \ (thrust - drag) / weight - print('spec energy rate:') - print('velocity, thrust, drag, weight') - print(velocity, thrust, drag, weight) def setup_partials(self): arange = np.arange(self.options['num_nodes']) diff --git a/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py b/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py index 825efac2c..87d2c6379 100644 --- a/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/compressibility_drag.py @@ -76,9 +76,6 @@ def compute(self, inputs, outputs): outputs['compress_drag_coeff'][idx_super] = cdc_super outputs['compress_drag_coeff'][idx_sub] = cdc_sub - print('compressibility drag:') - print('cdc_sub, cdc_super') - print(cdc_sub, cdc_super) def _compute_supersonic(self, inputs, outputs, idx): """ diff --git a/aviary/subsystems/aerodynamics/flops_based/drag.py b/aviary/subsystems/aerodynamics/flops_based/drag.py index 8cea0c931..0a1c4ee73 100644 --- a/aviary/subsystems/aerodynamics/flops_based/drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/drag.py @@ -66,9 +66,6 @@ def compute(self, inputs, outputs): q = inputs[Dynamic.Mission.DYNAMIC_PRESSURE] CD = inputs['drag_coefficient'] - print('simple drag') - print('M, q, CD') - print(M, q, CD) idx_sup = np.where(M >= 1.0) CD_scaled = CD * FCDSUB CD_scaled[idx_sup] = CD[idx_sup] * FCDSUP diff --git a/aviary/subsystems/aerodynamics/flops_based/induced_drag.py b/aviary/subsystems/aerodynamics/flops_based/induced_drag.py index e61278021..8ac3f4726 100644 --- a/aviary/subsystems/aerodynamics/flops_based/induced_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/induced_drag.py @@ -73,9 +73,6 @@ def compute(self, inputs, outputs): CL = 2.0 * lift / (Sref * gamma * P * mach ** 2) - print('induced drag:') - print('gamma, mach, lift, P, Sref, AR, span_efficiency_factor, SW25, TR') - print(gamma, mach, lift, P, Sref, AR, span_efficiency_factor, SW25, TR) redux, _ = aviary_options.get_item( Aircraft.Wing.SPAN_EFFICIENCY_REDUCTION, (False, None)) diff --git a/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py b/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py index ab251440e..adba49d48 100644 --- a/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/lift_dependent_drag.py @@ -115,9 +115,6 @@ def compute(self, inputs, outputs): gamma = self.options['gamma'] mach, lift, P, CLDES, MDES, Sref, AR, CAM, SW25, TC = inputs.values() - print('lift dependent drag:') - print('gamma, mach, lift, P, CLDES, MDES, Sref, AR, CAM, SW25, TC') - print(gamma, mach, lift, P, CLDES, MDES, Sref, AR, CAM, SW25, TC) FCDP = np.empty(nn, dtype=mach.dtype) DCDP = np.empty(nn, dtype=mach.dtype) dFCDP_dDELM = np.empty(nn, dtype=mach.dtype) diff --git a/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py b/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py index fb2ff3639..c7a9edd0f 100644 --- a/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py +++ b/aviary/subsystems/aerodynamics/flops_based/skin_friction_drag.py @@ -108,9 +108,6 @@ def compute(self, inputs, outputs): mission_wing_area = inputs[Aircraft.Wing.AREA] laminar_flow = np.any(lam_up > 0.0) or np.any(lam_low > 0.0) - print('skin friction:') - print('cf, Re, fineness, wetted_area, lam_up, lam_low, laminar_flow') - print(cf, Re, fineness, wetted_area, lam_up, lam_low, laminar_flow) if laminar_flow: laminar_upper = _calc_laminar_flow(lam_up) From 0e6f3fc573cf3e05661c44e2862a2f76d3d5ad1f Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Thu, 11 Jan 2024 08:35:22 -0800 Subject: [PATCH 12/34] Added SGM phases for detailed takeoff and landing --- aviary/mission/flops_based/ode/mission_ODE.py | 12 +- aviary/mission/flops_based/ode/range_rate.py | 1 - .../phases/time_integration_phases.py | 111 +++++++++++++++--- .../ode/time_integration_base_classes.py | 7 +- 4 files changed, 102 insertions(+), 29 deletions(-) diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index b48ad3dd5..67ba629dc 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -198,24 +198,21 @@ def setup(self): ], promotes_outputs=['initial_mass_residual']) - if analysis_scheme is AnalysisScheme.SHOOTING: + if analysis_scheme is AnalysisScheme.SHOOTING and False: 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 + 't_curr', Mission.Design.RESERVE_FUEL, Dynamic.Mission.MASS, Dynamic.Mission.RANGE, - Dynamic.Mission.THROTTLE, - 'required_lift', - 'load_factor', Dynamic.Mission.ALTITUDE, Dynamic.Mission.FLIGHT_PATH_ANGLE, ], input_units={ - 'required_lift': 'lbf', + 't_curr': 's', Dynamic.Mission.FLIGHT_PATH_ANGLE: 'deg', + Dynamic.Mission.RANGE: 'NM', }) self.add_subsystem( "dummy_comp", @@ -223,3 +220,4 @@ def setup(self): promotes_inputs=["*"],) self.set_input_defaults( Dynamic.Mission.RANGE, val=0, units='NM') + self.set_input_defaults('t_curr', val=0, units='s') diff --git a/aviary/mission/flops_based/ode/range_rate.py b/aviary/mission/flops_based/ode/range_rate.py index 7cf4b5693..9ef5d8fca 100644 --- a/aviary/mission/flops_based/ode/range_rate.py +++ b/aviary/mission/flops_based/ode/range_rate.py @@ -32,7 +32,6 @@ def compute(self, inputs, outputs): velocity = inputs[Dynamic.Mission.VELOCITY] climb_rate_2 = climb_rate**2 velocity_2 = velocity**2 - print(climb_rate, velocity) if (climb_rate_2 >= velocity_2).any(): raise om.AnalysisError( "WARNING: climb rate exceeds velocity (range_rate.py)") diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index 8cd771e22..05c8651f9 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -3,6 +3,8 @@ import openmdao.api as om 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.flops_based.phases.cruise_phase import Cruise from aviary.mission.gasp_based.phases.time_integration_traj import FlexibleTraj from aviary.mission.gasp_based.ode.time_integration_base_classes import SimuPyProblem @@ -16,12 +18,17 @@ class SGMHeightEnergy(SimuPyProblem): - def __init__(self, - ode_args={}, - simupy_args={},): + def __init__( + self, + phase_name='cruise', + distance_trigger_units='NM', + ode_args={}, + simupy_args={}, + ): super().__init__(MissionODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), + output_names=[], state_names=[ Dynamic.Mission.MASS, Dynamic.Mission.RANGE, @@ -31,6 +38,8 @@ def __init__(self, Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args) + self.phase_name = phase_name + self.distance_trigger_units = distance_trigger_units self.event_channel_names = [ # Dynamic.Mission.DISTANCE, Dynamic.Mission.MASS, @@ -39,19 +48,89 @@ def __init__(self, def event_equation_function(self, t, x): self.output_equation_function(t, x) - 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() + # distance = self.get_val(Dynamic.Mission.RANGE, + # units=self.distance_trigger_units).squeeze() + # distance_trigger = self.get_val( + # "distance_trigger", units=self.distance_trigger_units).squeeze() current_mass = self.get_val(Dynamic.Mission.MASS, units="lbm").squeeze() mass_trigger = 150000 - return np.array([ current_mass - mass_trigger ]) +class SGMDetailedTakeoff(SimuPyProblem): + def __init__( + self, + phase_name='detailed_takeoff', + ode_args={}, + simupy_args={}, + ): + super().__init__(TakeoffODE( + analysis_scheme=AnalysisScheme.SHOOTING, + **ode_args), + output_names=[], + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.RANGE, + Dynamic.Mission.ALTITUDE, + ], + alternate_state_rate_names={ + Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + **simupy_args) + + self.phase_name = phase_name + self.event_channel_names = [ + Dynamic.Mission.ALTITUDE, + ] + self.num_events = len(self.event_channel_names) + + def event_equation_function(self, t, x): + self.output_equation_function(t, x) + current_alt = self.get_val(Dynamic.Mission.ALTITUDE, units="ft").squeeze() + alt_trigger = 50 + return np.array([ + current_alt - alt_trigger + # maybe mach + ]) + + +class SGMDetailedLanding(SimuPyProblem): + def __init__( + self, + phase_name='detailed_landing', + ode_args={}, + simupy_args={}, + ): + super().__init__(LandingODE( + analysis_scheme=AnalysisScheme.SHOOTING, + **ode_args), + output_names=[], + state_names=[ + Dynamic.Mission.MASS, + Dynamic.Mission.RANGE, + Dynamic.Mission.ALTITUDE, + ], + alternate_state_rate_names={ + Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + **simupy_args) + + self.phase_name = phase_name + self.event_channel_names = [ + Dynamic.Mission.ALTITUDE, + ] + self.num_events = len(self.event_channel_names) + + def event_equation_function(self, t, x): + self.output_equation_function(t, x) + current_alt = self.get_val(Dynamic.Mission.ALTITUDE, units="ft").squeeze() + alt_trigger = 0 + return np.array([ + current_alt - alt_trigger + ]) + + def test_phase(phases, ode_args_tab): prob = om.Problem() prob.driver = om.pyOptSparseDriver() @@ -91,13 +170,15 @@ def test_phase(phases, ode_args_tab): promotes_inputs=['aircraft:*', 'mission:*'], promotes_outputs=['aircraft:*', 'mission:*'] ) - prob.model.add_subsystem('traj', traj,) # promotes=['aircraft:*','mission:*'] - prob.model.promotes('traj', inputs=[ - Aircraft.Wing.CHARACTERISTIC_LENGTH, - Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, - Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, - Aircraft.Fuselage.CHARACTERISTIC_LENGTH, - Aircraft.Nacelle.CHARACTERISTIC_LENGTH,]) + prob.model.add_subsystem('traj', traj, + promotes=['aircraft:*', 'mission:*'] + ) + # prob.model.promotes('traj', inputs=[ + # Aircraft.Wing.CHARACTERISTIC_LENGTH, + # Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, + # Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, + # Aircraft.Fuselage.CHARACTERISTIC_LENGTH, + # Aircraft.Nacelle.CHARACTERISTIC_LENGTH,]) prob.model.add_subsystem( "fuel_obj", 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 0a923ae40..59e11f10e 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -179,7 +179,7 @@ def __init__( for output_name in output_names ] - if include_state_outputs: + if include_state_outputs or output_names == []: # prevent empty outputs output_names = state_names + output_names self.output_units = self.state_units + self.output_units @@ -227,8 +227,6 @@ def state(self): @state.setter def state(self, value): - print('state.setter: self.state, value') - print(self.state, value) if np.all(self.state == value): return for state_name, elem_val, unit in zip( @@ -414,9 +412,6 @@ def setup_params( if '.' in prom_name: continue traj_promote_initial_input[prom_name] = data - # print(traj_promote_initial_input) - # print(len(traj_promote_initial_input)) - # exit() self.traj_promote_initial_input = { **self.options["param_dict"], **traj_promote_initial_input} From 7349550443b7d2f6297bfa1a7c163f9cd604defa Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Wed, 17 Jan 2024 09:44:33 -0800 Subject: [PATCH 13/34] Switched from mission_ODE to simple_mission_ODE --- .../flops_based/ode/simple_mission_ODE.py | 40 +++++++++++++++++++ .../phases/time_integration_phases.py | 4 +- 2 files changed, 42 insertions(+), 2 deletions(-) diff --git a/aviary/mission/flops_based/ode/simple_mission_ODE.py b/aviary/mission/flops_based/ode/simple_mission_ODE.py index 58554956b..0df9dd42f 100644 --- a/aviary/mission/flops_based/ode/simple_mission_ODE.py +++ b/aviary/mission/flops_based/ode/simple_mission_ODE.py @@ -2,12 +2,14 @@ import openmdao.api as om from dymos.models.atmosphere import USatm1976Comp +from aviary.mission.gasp_based.ode.time_integration_base_classes import add_SGM_required_inputs from aviary.mission.flops_based.ode.simple_mission_EOM import MissionEOM 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_in import VariablesIn +from aviary.variable_info.enums import AnalysisScheme class ExternalSubsystemGroup(om.Group): @@ -43,15 +45,29 @@ def initialize(self): 'throttle_enforcement', default='path', 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.RANGE: {'units': 'm'}, + } + add_SGM_required_inputs(self, SGM_required_inputs) + self.add_subsystem( 'input_port', VariablesIn(aviary_options=aviary_options, @@ -204,3 +220,27 @@ def setup(self): self.linear_solver = om.DirectSolver(assemble_jac=True) self.nonlinear_solver.options['err_on_non_converge'] = True self.nonlinear_solver.options['iprint'] = 2 + + if analysis_scheme is AnalysisScheme.SHOOTING and False: + from aviary.utils.functions import create_printcomp + dummy_comp = create_printcomp( + all_inputs=[ + 't_curr', + Mission.Design.RESERVE_FUEL, + Dynamic.Mission.MASS, + Dynamic.Mission.RANGE, + Dynamic.Mission.ALTITUDE, + Dynamic.Mission.FLIGHT_PATH_ANGLE, + ], + input_units={ + 't_curr': 's', + Dynamic.Mission.FLIGHT_PATH_ANGLE: 'deg', + Dynamic.Mission.RANGE: 'NM', + }) + self.add_subsystem( + "dummy_comp", + dummy_comp(), + promotes_inputs=["*"],) + self.set_input_defaults( + Dynamic.Mission.RANGE, val=0, units='NM') + self.set_input_defaults('t_curr', val=0, units='s') diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index 05c8651f9..a879bb3e5 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -2,7 +2,7 @@ import dymos as dm import openmdao.api as om -from aviary.mission.flops_based.ode.mission_ODE import MissionODE +from aviary.mission.flops_based.ode.simple_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.flops_based.phases.cruise_phase import Cruise @@ -233,7 +233,7 @@ def test_phase(phases, ode_args_tab): if __name__ == '__main__': - from aviary.interface.default_phase_info.flops import phase_info, aero, prop, geom + from aviary.interface.default_phase_info.simple import phase_info, aero, prop, geom from aviary.subsystems.propulsion.engine_deck import EngineDeck from aviary.variable_info.variables import Aircraft, Mission, Dynamic from aviary.utils.process_input_decks import create_vehicle From f4c7017a16964b3e8c708e146b5aca7f34f5915a Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Thu, 18 Jan 2024 15:38:52 -0800 Subject: [PATCH 14/34] Converting simple_mission to SGM Also added some error checking to automatic state detections --- aviary/mission/flops_based/ode/mission_EOM.py | 6 ++-- .../flops_based/ode/simple_mission_ODE.py | 35 ++++++++++++++++++- .../phases/time_integration_phases.py | 10 ++++-- .../ode/time_integration_base_classes.py | 14 ++++++-- 4 files changed, 57 insertions(+), 8 deletions(-) diff --git a/aviary/mission/flops_based/ode/mission_EOM.py b/aviary/mission/flops_based/ode/mission_EOM.py index c679ed093..97ba87013 100644 --- a/aviary/mission/flops_based/ode/mission_EOM.py +++ b/aviary/mission/flops_based/ode/mission_EOM.py @@ -4,15 +4,12 @@ from aviary.mission.flops_based.ode.range_rate import RangeRate from aviary.mission.flops_based.ode.specific_energy_rate import SpecificEnergyRate from aviary.variable_info.variables import Dynamic -from aviary.variable_info.enums import AnalysisScheme class MissionEOM(om.Group): def initialize(self): self.options.declare('num_nodes', types=int, desc='Number of nodes to be evaluated in the RHS') - self.options.declare("analysis_scheme", types=AnalysisScheme, default=AnalysisScheme.COLLOCATION, - desc="The analysis method that will be used to close the trajectory; for example collocation or time integration") def setup(self): nn = self.options['num_nodes'] @@ -22,16 +19,19 @@ def setup(self): promotes_inputs=[Dynamic.Mission.VELOCITY, Dynamic.Mission.MASS, Dynamic.Mission.THRUST_TOTAL, Dynamic.Mission.DRAG], promotes_outputs=[Dynamic.Mission.SPECIFIC_ENERGY_RATE]) + self.add_subsystem(name=Dynamic.Mission.ALTITUDE_RATE, subsys=AltitudeRate(num_nodes=nn), promotes_inputs=[Dynamic.Mission.SPECIFIC_ENERGY_RATE, Dynamic.Mission.VELOCITY_RATE, Dynamic.Mission.VELOCITY], promotes_outputs=[Dynamic.Mission.ALTITUDE_RATE]) + self.add_subsystem(name='groundspeed', subsys=RangeRate(num_nodes=nn), promotes_inputs=[ Dynamic.Mission.ALTITUDE_RATE, Dynamic.Mission.VELOCITY], promotes_outputs=[Dynamic.Mission.RANGE_RATE]) + self.add_subsystem(name='excess_specific_power', subsys=SpecificEnergyRate(num_nodes=nn), promotes_inputs=[(Dynamic.Mission.THRUST_TOTAL, Dynamic.Mission.THRUST_MAX_TOTAL), diff --git a/aviary/mission/flops_based/ode/simple_mission_ODE.py b/aviary/mission/flops_based/ode/simple_mission_ODE.py index 0df9dd42f..8adca7d5b 100644 --- a/aviary/mission/flops_based/ode/simple_mission_ODE.py +++ b/aviary/mission/flops_based/ode/simple_mission_ODE.py @@ -7,7 +7,7 @@ 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 @@ -183,9 +183,42 @@ def setup(self): promotes_inputs=['*'], promotes_outputs=['*']) + self.set_input_defaults(Mission.Design.GROSS_MASS, val=1, units='kg') + self.set_input_defaults( + Aircraft.Fuselage.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.Fuselage.FINENESS, val=1, units='unitless') + self.set_input_defaults(Aircraft.Fuselage.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults( + Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.VerticalTail.FINENESS, val=1, units='unitless') + self.set_input_defaults(Aircraft.VerticalTail.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults( + Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.HorizontalTail.FINENESS, + val=1, units='unitless') + self.set_input_defaults( + Aircraft.HorizontalTail.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults(Aircraft.Wing.CHARACTERISTIC_LENGTH, val=1, units='ft') + self.set_input_defaults(Aircraft.Wing.FINENESS, val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.WETTED_AREA, val=1, units='ft**2') + self.set_input_defaults( + Aircraft.Wing.SPAN_EFFICIENCY_FACTOR, val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.TAPER_RATIO, val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.THICKNESS_TO_CHORD, + val=1, units='unitless') + self.set_input_defaults(Aircraft.Wing.SWEEP, val=1, units='deg') + self.set_input_defaults(Aircraft.Wing.ASPECT_RATIO, val=1, units='unitless') + self.set_input_defaults( + Aircraft.Design.LIFT_DEPENDENT_DRAG_COEFF_FACTOR, val=1, units='unitless') + self.set_input_defaults( + Aircraft.Design.ZERO_LIFT_DRAG_COEFF_FACTOR, val=1, units='unitless') + + 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') diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index a879bb3e5..28afbdff4 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -33,6 +33,11 @@ def __init__( Dynamic.Mission.MASS, Dynamic.Mission.RANGE, Dynamic.Mission.ALTITUDE, + ], + state_rate_units=[ + 'lbm/s', + 'm/s', + 'm/s', ], alternate_state_rate_names={ Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, @@ -214,8 +219,9 @@ def test_phase(phases, ode_args_tab): prob.set_val("traj.altitude_initial", val=35000, units="ft") prob.set_val("traj.mass_initial", val=171000, units="lbm") prob.set_val("traj.range_initial", val=0, units="NM") - prob.set_val("traj.velocity", val=472, units="kn") - prob.set_val("traj.velocity_rate", val=0, units="m/s**2") + prob.set_val("traj.mach", val=.8, units="unitless") + # prob.set_val("traj.velocity", val=472, units="kn") + # prob.set_val("traj.velocity_rate", val=0, units="m/s**2") # try: prob.run_model() 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 59e11f10e..f67f1d57c 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -26,11 +26,11 @@ def __init__( time_independent=False, t_name="t_curr", state_names=None, - alternate_state_names={}, + alternate_state_names=None, blocked_state_names=None, state_units=None, state_rate_names=None, - alternate_state_rate_names={}, + alternate_state_rate_names=None, state_rate_units=None, parameter_names=None, parameter_units=None, @@ -95,6 +95,10 @@ def __init__( outputs = [val['prom_name'] for key, val in data] data = prob.model.list_inputs(prom_name=True, val=False, out_stream=None) inputs = [val['prom_name'] for key, val in data] + if alternate_state_names is None: + alternate_state_names = {} + if alternate_state_rate_names is None: + alternate_state_rate_names = {} if state_names is None: state_names = [ @@ -133,6 +137,9 @@ def __init__( units=True, ) ] + if len(self.state_rate_units) != len(state_rate_names): + raise RuntimeError( + 'state_rate_units could not be auto generated, declare them explicitly.') if state_units is not None: self.state_units = state_units @@ -141,6 +148,9 @@ def __init__( units.simplify_unit("s*" + rate_unit) for rate_unit in self.state_rate_units ] + if len(self.state_units) != len(state_names): + raise RuntimeError( + 'state_units could not be auto generated, declare them explicitly.') if parameter_names is None: parameter_names = [ From 218e803ae8e9bcda4a817a1707c3b409ba564eb6 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 23 Jan 2024 16:27:03 -0800 Subject: [PATCH 15/34] prepared for testing detailed takeoff and landing --- aviary/mission/flops_based/ode/landing_ode.py | 2 +- aviary/mission/flops_based/ode/mission_ODE.py | 23 +- .../flops_based/ode/simple_mission_EOM.py | 46 --- .../flops_based/ode/simple_mission_ODE.py | 279 ------------------ aviary/mission/flops_based/ode/takeoff_ode.py | 2 +- .../phases/time_integration_phases.py | 96 +++--- .../ode/time_integration_base_classes.py | 23 +- 7 files changed, 88 insertions(+), 383 deletions(-) delete mode 100644 aviary/mission/flops_based/ode/simple_mission_EOM.py delete mode 100644 aviary/mission/flops_based/ode/simple_mission_ODE.py diff --git a/aviary/mission/flops_based/ode/landing_ode.py b/aviary/mission/flops_based/ode/landing_ode.py index 734623ae1..48827e681 100644 --- a/aviary/mission/flops_based/ode/landing_ode.py +++ b/aviary/mission/flops_based/ode/landing_ode.py @@ -81,7 +81,7 @@ def setup(self): if analysis_scheme is AnalysisScheme.SHOOTING: SGM_required_inputs = { 't_curr': {'units': 's'}, - Dynamic.RANGE: {'units': 'm'}, + Dynamic.Mission.DISTANCE: {'units': 'm'}, } add_SGM_required_inputs(self, SGM_required_inputs) diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index 2257d0b5e..600aa83f0 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -3,7 +3,7 @@ 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 +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 @@ -65,7 +65,7 @@ def setup(self): if analysis_scheme is AnalysisScheme.SHOOTING: SGM_required_inputs = { 't_curr': {'units': 's'}, - Dynamic.Mission.RANGE: {'units': 'm'}, + Dynamic.Mission.DISTANCE: {'units': 'm'}, } add_SGM_required_inputs(self, SGM_required_inputs) @@ -214,9 +214,12 @@ def setup(self): self.set_input_defaults( Aircraft.Design.ZERO_LIFT_DRAG_COEFF_FACTOR, val=1, units='unitless') + 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') @@ -243,6 +246,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 + if analysis_scheme is AnalysisScheme.SHOOTING and False: from aviary.utils.functions import create_printcomp dummy_comp = create_printcomp( @@ -250,21 +261,21 @@ def setup(self): 't_curr', Mission.Design.RESERVE_FUEL, Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, + Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, Dynamic.Mission.FLIGHT_PATH_ANGLE, ], input_units={ 't_curr': 's', Dynamic.Mission.FLIGHT_PATH_ANGLE: 'deg', - Dynamic.Mission.RANGE: 'NM', + Dynamic.Mission.DISTANCE: 'NM', }) self.add_subsystem( "dummy_comp", dummy_comp(), promotes_inputs=["*"],) self.set_input_defaults( - Dynamic.Mission.RANGE, val=0, units='NM') + Dynamic.Mission.DISTANCE, val=0, units='NM') self.set_input_defaults('t_curr', val=0, units='s') self.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, @@ -274,4 +285,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/simple_mission_EOM.py b/aviary/mission/flops_based/ode/simple_mission_EOM.py deleted file mode 100644 index 55bc58889..000000000 --- a/aviary/mission/flops_based/ode/simple_mission_EOM.py +++ /dev/null @@ -1,46 +0,0 @@ -import openmdao.api as om - -from aviary.variable_info.variables import Dynamic -from aviary.mission.flops_based.ode.altitude_rate import AltitudeRate -from aviary.mission.flops_based.ode.range_rate import RangeRate -from aviary.mission.flops_based.ode.specific_energy_rate import SpecificEnergyRate -from aviary.mission.flops_based.ode.required_thrust import RequiredThrust - - -class MissionEOM(om.Group): - def initialize(self): - self.options.declare('num_nodes', types=int, - desc='Number of nodes to be evaluated in the RHS') - - def setup(self): - nn = self.options['num_nodes'] - - self.add_subsystem(name='required_thrust', - subsys=RequiredThrust(num_nodes=nn), - promotes_inputs=[Dynamic.Mission.DRAG, Dynamic.Mission.ALTITUDE_RATE, - Dynamic.Mission.VELOCITY, Dynamic.Mission.VELOCITY_RATE, Dynamic.Mission.MASS], - promotes_outputs=['thrust_required']) - - self.add_subsystem(name='groundspeed', - subsys=RangeRate(num_nodes=nn), - promotes_inputs=[ - Dynamic.Mission.ALTITUDE_RATE, Dynamic.Mission.VELOCITY], - promotes_outputs=[Dynamic.Mission.DISTANCE_RATE]) - - self.add_subsystem(name='excess_specific_power', - subsys=SpecificEnergyRate(num_nodes=nn), - promotes_inputs=[(Dynamic.Mission.THRUST_TOTAL, Dynamic.Mission.THRUST_MAX_TOTAL), - Dynamic.Mission.VELOCITY, Dynamic.Mission.MASS, Dynamic.Mission.DRAG], - promotes_outputs=[(Dynamic.Mission.SPECIFIC_ENERGY_RATE, Dynamic.Mission.SPECIFIC_ENERGY_RATE_EXCESS)]) - self.add_subsystem( - name=Dynamic.Mission.ALTITUDE_RATE_MAX, - subsys=AltitudeRate( - num_nodes=nn), - promotes_inputs=[ - (Dynamic.Mission.SPECIFIC_ENERGY_RATE, - Dynamic.Mission.SPECIFIC_ENERGY_RATE_EXCESS), - Dynamic.Mission.VELOCITY_RATE, - Dynamic.Mission.VELOCITY], - promotes_outputs=[ - (Dynamic.Mission.ALTITUDE_RATE, - Dynamic.Mission.ALTITUDE_RATE_MAX)]) diff --git a/aviary/mission/flops_based/ode/simple_mission_ODE.py b/aviary/mission/flops_based/ode/simple_mission_ODE.py deleted file mode 100644 index dc224eb43..000000000 --- a/aviary/mission/flops_based/ode/simple_mission_ODE.py +++ /dev/null @@ -1,279 +0,0 @@ -import numpy as np -import openmdao.api as om -from dymos.models.atmosphere import USatm1976Comp - -from aviary.mission.gasp_based.ode.time_integration_base_classes import add_SGM_required_inputs -from aviary.mission.flops_based.ode.simple_mission_EOM import MissionEOM -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, Aircraft -from aviary.variable_info.variables_in import VariablesIn -from aviary.variable_info.enums import AnalysisScheme - - -class ExternalSubsystemGroup(om.Group): - def configure(self): - promote_aircraft_and_mission_vars(self) - - -class MissionODE(om.Group): - def initialize(self): - self.options.declare( - 'num_nodes', types=int, - desc='Number of nodes to be evaluated in the RHS') - self.options.declare( - 'subsystem_options', types=dict, default={}, - desc='dictionary of parameters to be passed to the subsystem builders') - self.options.declare( - 'aviary_options', types=AviaryValues, - desc='collection of Aircraft/Mission specific options') - self.options.declare( - 'core_subsystems', - desc='list of core subsystem builder instances to be added to the ODE' - ) - self.options.declare( - 'external_subsystems', default=[], - desc='list of external subsystem builder instances to be added to the ODE') - self.options.declare( - 'meta_data', default=_MetaData, - desc='metadata associated with the variables to be passed into the ODE') - self.options.declare( - 'use_actual_takeoff_mass', default=False, - desc='flag to use actual takeoff mass in the climb phase, otherwise assume 100 kg fuel burn') - self.options.declare( - 'throttle_enforcement', default='path', - 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.RANGE: {'units': 'm'}, - } - add_SGM_required_inputs(self, SGM_required_inputs) - - self.add_subsystem( - 'input_port', - VariablesIn(aviary_options=aviary_options, - meta_data=self.options['meta_data'], - context='mission'), - promotes_inputs=['*'], - promotes_outputs=['*']) - self.add_subsystem( - name='atmosphere', - subsys=USatm1976Comp(num_nodes=nn), - promotes_inputs=[('h', Dynamic.Mission.ALTITUDE)], - promotes_outputs=[ - ('sos', Dynamic.Mission.SPEED_OF_SOUND), ('rho', Dynamic.Mission.DENSITY), - ('temp', Dynamic.Mission.TEMPERATURE), ('pres', Dynamic.Mission.STATIC_PRESSURE)]) - - # add an execcomp to compute velocity based off mach and sos - self.add_subsystem( - name='velocity_comp', - subsys=om.ExecComp( - 'velocity = mach * sos', - mach={'units': 'unitless', 'shape': (nn,)}, - sos={'units': 'm/s', 'shape': (nn,)}, - velocity={'units': 'm/s', 'shape': (nn,)}, - has_diag_partials=True, - ), - promotes_inputs=[('mach', Dynamic.Mission.MACH), - ('sos', Dynamic.Mission.SPEED_OF_SOUND)], - promotes_outputs=[('velocity', Dynamic.Mission.VELOCITY)], - ) - - # add execcomp to compute velocity_rate based off mach_rate and sos - self.add_subsystem( - name='velocity_rate_comp', - subsys=om.ExecComp( - 'velocity_rate = mach_rate * sos', - mach_rate={'units': 'unitless/s', 'shape': (nn,)}, - sos={'units': 'm/s', 'shape': (nn,)}, - velocity_rate={'units': 'm/s**2', 'shape': (nn,)}, - has_diag_partials=True, - ), - promotes_inputs=[('mach_rate', Dynamic.Mission.MACH_RATE), - ('sos', Dynamic.Mission.SPEED_OF_SOUND)], - promotes_outputs=[('velocity_rate', Dynamic.Mission.VELOCITY_RATE)], - ) - - base_options = {'num_nodes': nn, 'aviary_inputs': aviary_options} - - for subsystem in core_subsystems: - # check if subsystem_options has entry for a subsystem of this name - if subsystem.name in subsystem_options: - kwargs = subsystem_options[subsystem.name] - else: - kwargs = {} - - kwargs.update(base_options) - system = subsystem.build_mission(**kwargs) - - if system is not None: - self.add_subsystem(subsystem.name, - system, - promotes_inputs=subsystem.mission_inputs(**kwargs), - promotes_outputs=subsystem.mission_outputs(**kwargs)) - - # Create a lightly modified version of an OM group to add external subsystems - # to the ODE with a special configure() method that promotes - # all aircraft:* and mission:* variables to the ODE. - external_subsystem_group = ExternalSubsystemGroup() - add_subsystem_group = False - - for subsystem in self.options['external_subsystems']: - subsystem_mission = subsystem.build_mission( - num_nodes=nn, aviary_inputs=aviary_options) - if subsystem_mission is not None: - add_subsystem_group = True - external_subsystem_group.add_subsystem(subsystem.name, subsystem_mission) - - # Only add the external subsystem group if it has at least one subsystem. - # Without this logic there'd be an empty OM group added to the ODE. - if add_subsystem_group: - self.add_subsystem( - name='external_subsystems', - subsys=external_subsystem_group, - promotes_inputs=['*'], - promotes_outputs=['*']) - - self.add_subsystem( - name='mission_EOM', - subsys=MissionEOM(num_nodes=nn), - promotes_inputs=[ - Dynamic.Mission.VELOCITY, Dynamic.Mission.MASS, - Dynamic.Mission.THRUST_MAX_TOTAL, - Dynamic.Mission.DRAG, Dynamic.Mission.ALTITUDE_RATE, - Dynamic.Mission.VELOCITY_RATE], - promotes_outputs=[ - Dynamic.Mission.SPECIFIC_ENERGY_RATE_EXCESS, - Dynamic.Mission.ALTITUDE_RATE_MAX, - Dynamic.Mission.DISTANCE_RATE, - 'thrust_required', - ]) - - # add a balance comp to compute throttle based on the altitude rate - self.add_subsystem(name='throttle_balance', - subsys=om.BalanceComp(name=Dynamic.Mission.THROTTLE, - units="unitless", - val=np.ones(nn), - lhs_name='thrust_required', - rhs_name=Dynamic.Mission.THRUST_TOTAL, - eq_units="lbf", - normalize=False, - lower=0.0 if options['throttle_enforcement'] == 'bounded' else None, - upper=1.0 if options['throttle_enforcement'] == 'bounded' else None, - ), - promotes_inputs=['*'], - promotes_outputs=['*']) - - self.set_input_defaults(Mission.Design.GROSS_MASS, val=1, units='kg') - self.set_input_defaults( - Aircraft.Fuselage.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.Fuselage.FINENESS, val=1, units='unitless') - self.set_input_defaults(Aircraft.Fuselage.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults( - Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.VerticalTail.FINENESS, val=1, units='unitless') - self.set_input_defaults(Aircraft.VerticalTail.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults( - Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.HorizontalTail.FINENESS, - val=1, units='unitless') - self.set_input_defaults( - Aircraft.HorizontalTail.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults(Aircraft.Wing.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.Wing.FINENESS, val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults( - Aircraft.Wing.SPAN_EFFICIENCY_FACTOR, val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.TAPER_RATIO, val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.THICKNESS_TO_CHORD, - val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.SWEEP, val=1, units='deg') - self.set_input_defaults(Aircraft.Wing.ASPECT_RATIO, val=1, units='unitless') - self.set_input_defaults( - Aircraft.Design.LIFT_DEPENDENT_DRAG_COEFF_FACTOR, val=1, units='unitless') - self.set_input_defaults( - Aircraft.Design.ZERO_LIFT_DRAG_COEFF_FACTOR, val=1, units='unitless') - - 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') - - if options['use_actual_takeoff_mass']: - exec_comp_string = 'initial_mass_residual = initial_mass - mass[0]' - initial_mass_string = Mission.Takeoff.FINAL_MASS - else: - exec_comp_string = 'initial_mass_residual = initial_mass - mass[0] - 100.' - initial_mass_string = Mission.Design.GROSS_MASS - - # Experimental: Add a component to constrain the initial mass to be equal to design gross weight. - initial_mass_residual_constraint = om.ExecComp( - exec_comp_string, - initial_mass={'units': 'kg'}, - mass={'units': 'kg', 'shape': (nn,)}, - initial_mass_residual={'units': 'kg'}, - ) - - self.add_subsystem('initial_mass_residual_constraint', initial_mass_residual_constraint, - promotes_inputs=[ - ('initial_mass', initial_mass_string), - ('mass', Dynamic.Mission.MASS) - ], - promotes_outputs=['initial_mass_residual']) - - self.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, - atol=1.0e-10, - rtol=1.0e-10, - ) - 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 - - if analysis_scheme is AnalysisScheme.SHOOTING and False: - from aviary.utils.functions import create_printcomp - dummy_comp = create_printcomp( - all_inputs=[ - 't_curr', - Mission.Design.RESERVE_FUEL, - Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, - Dynamic.Mission.ALTITUDE, - Dynamic.Mission.FLIGHT_PATH_ANGLE, - ], - input_units={ - 't_curr': 's', - Dynamic.Mission.FLIGHT_PATH_ANGLE: 'deg', - Dynamic.Mission.RANGE: 'NM', - }) - self.add_subsystem( - "dummy_comp", - dummy_comp(), - promotes_inputs=["*"],) - self.set_input_defaults( - Dynamic.Mission.RANGE, val=0, units='NM') - self.set_input_defaults('t_curr', val=0, units='s') diff --git a/aviary/mission/flops_based/ode/takeoff_ode.py b/aviary/mission/flops_based/ode/takeoff_ode.py index bfb55e790..e7292603b 100644 --- a/aviary/mission/flops_based/ode/takeoff_ode.py +++ b/aviary/mission/flops_based/ode/takeoff_ode.py @@ -79,7 +79,7 @@ def setup(self): if analysis_scheme is AnalysisScheme.SHOOTING: SGM_required_inputs = { 't_curr': {'units': 's'}, - Dynamic.RANGE: {'units': 'm'}, + Dynamic.Mission.DISTANCE: {'units': 'm'}, } add_SGM_required_inputs(self, SGM_required_inputs) diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index 28afbdff4..dc098ca88 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -1,14 +1,11 @@ import numpy as np -import dymos as dm import openmdao.api as om -from aviary.mission.flops_based.ode.simple_mission_ODE import MissionODE +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.flops_based.phases.cruise_phase import Cruise from aviary.mission.gasp_based.phases.time_integration_traj import FlexibleTraj 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 from aviary.variable_info.variables import Aircraft, Dynamic, Mission from aviary.variable_info.variables_in import VariablesIn @@ -28,19 +25,19 @@ def __init__( super().__init__(MissionODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), - output_names=[], - state_names=[ + outputs=[], + states=[ Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, + Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, ], - state_rate_units=[ - 'lbm/s', - 'm/s', - 'm/s', - ], + # state_rate_units=[ + # 'lbm/s', + # 'm/s', + # 'm/s', + # ], alternate_state_rate_names={ - Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + Dynamic.Mission.MASS: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args) self.phase_name = phase_name @@ -53,7 +50,7 @@ def __init__( def event_equation_function(self, t, x): self.output_equation_function(t, x) - # distance = self.get_val(Dynamic.Mission.RANGE, + # 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() @@ -75,14 +72,14 @@ def __init__( super().__init__(TakeoffODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), - output_names=[], - state_names=[ + outputs=[], + states=[ Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, + Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, ], alternate_state_rate_names={ - Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + Dynamic.Mission.MASS: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args) self.phase_name = phase_name @@ -111,14 +108,14 @@ def __init__( super().__init__(LandingODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), - output_names=[], - state_names=[ + outputs=[], + states=[ Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, + Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, ], alternate_state_rate_names={ - Dynamic.Mission.MASS_RATE: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, + Dynamic.Mission.MASS: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args) self.phase_name = phase_name @@ -151,19 +148,12 @@ def test_phase(phases, ode_args_tab): traj = FlexibleTraj( Phases=phases, promote_all_auto_ivc=True, - # traj_promote_initial_input={ - # Aircraft.Wing.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - # Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - # Aircraft.VerticalTail.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - # Aircraft.Fuselage.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - # Aircraft.Nacelle.CHARACTERISTIC_LENGTH: {'units': 'ft'}, - # }, traj_final_state_output=[Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, + Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE], traj_initial_state_input=[ Dynamic.Mission.MASS, - Dynamic.Mission.RANGE, + Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, ], ) @@ -178,12 +168,6 @@ def test_phase(phases, ode_args_tab): prob.model.add_subsystem('traj', traj, promotes=['aircraft:*', 'mission:*'] ) - # prob.model.promotes('traj', inputs=[ - # Aircraft.Wing.CHARACTERISTIC_LENGTH, - # Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, - # Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, - # Aircraft.Fuselage.CHARACTERISTIC_LENGTH, - # Aircraft.Nacelle.CHARACTERISTIC_LENGTH,]) prob.model.add_subsystem( "fuel_obj", @@ -216,10 +200,28 @@ def test_phase(phases, ode_args_tab): warnings.simplefilter("ignore", om.PromotionWarning) prob.setup() - prob.set_val("traj.altitude_initial", val=35000, units="ft") - prob.set_val("traj.mass_initial", val=171000, units="lbm") - prob.set_val("traj.range_initial", val=0, units="NM") - prob.set_val("traj.mach", val=.8, units="unitless") + + 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': 0, 'units': "m/s"}, + } + 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"}, + } + 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"}, + } + + for key, val in initial_values_cruise.items(): + prob.set_val(key, **val) # prob.set_val("traj.velocity", val=472, units="kn") # prob.set_val("traj.velocity_rate", val=0, units="m/s**2") @@ -233,13 +235,13 @@ def test_phase(phases, ode_args_tab): with open('output_list.txt', 'w') as outfile: prob.model.list_outputs(out_stream=outfile) - final_range = prob.get_val('traj.range_final', units='NM')[0] + final_distance = prob.get_val('traj.distance_final', units='NM')[0] final_mass = prob.get_val('traj.mass_final', units='lbm')[0] - print(final_range, final_mass) + print(final_distance, final_mass) if __name__ == '__main__': - from aviary.interface.default_phase_info.simple import phase_info, aero, prop, geom + from aviary.interface.default_phase_info.height_energy import phase_info, aero, prop, geom from aviary.subsystems.propulsion.engine_deck import EngineDeck from aviary.variable_info.variables import Aircraft, Mission, Dynamic from aviary.utils.process_input_decks import create_vehicle @@ -253,12 +255,18 @@ def test_phase(phases, ode_args_tab): aviary_inputs.set_val('debug_mode', False) 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_tab = dict(aviary_options=aviary_inputs, core_subsystems=core_subsystems) engine = EngineDeck(options=aviary_inputs) preprocess_propulsion(aviary_inputs, [engine]) ode_args_tab['num_nodes'] = 1 ode_args_tab['subsystem_options'] = {'core_aerodynamics': {'method': 'computed'}} + # ode_args_tab['friction_key'] = Mission.Takeoff.ROLLING_FRICTION_COEFFICIENT + # ode_args_tab['friction_key'] = Mission.Takeoff.BRAKING_FRICTION_COEFFICIENT phase_kwargs = dict( ode_args=ode_args_tab, simupy_args=dict( 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 b3a8d0e17..6b3a59afc 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -9,16 +9,27 @@ 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()), ) +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()), + ) + # Subproblem used as a basis for forward in time integration phases. + + class SimuPyProblem(SimulationMixin): def __init__( self, @@ -159,7 +170,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 } @@ -171,7 +182,7 @@ 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 @@ -186,7 +197,7 @@ def __init__( # TODO: add defensive checks to make sure dimensions match in both setup and # calls - if DEBUG or True: + if DEBUG: om.n2(prob, outfile="n2_simupy_problem.html", show_browser=False) with open('input_list_simupy.txt', 'w') as outfile: prob.model.list_inputs(out_stream=outfile,) From f2fb553f7c53f4c8cdbd4f67fd6971fc810484b1 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Wed, 31 Jan 2024 15:43:46 -0800 Subject: [PATCH 16/34] cleaned up comments and removed unused values also added example of customizing an event function --- .../phases/time_integration_phases.py | 33 ++++++++----------- aviary/mission/ode/specific_energy_rate.py | 3 -- 2 files changed, 14 insertions(+), 22 deletions(-) diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index dc098ca88..d7a9eb9f0 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -18,7 +18,6 @@ class SGMHeightEnergy(SimuPyProblem): def __init__( self, phase_name='cruise', - distance_trigger_units='NM', ode_args={}, simupy_args={}, ): @@ -31,29 +30,18 @@ def __init__( Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, ], - # state_rate_units=[ - # 'lbm/s', - # 'm/s', - # 'm/s', - # ], alternate_state_rate_names={ Dynamic.Mission.MASS: Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL}, **simupy_args) self.phase_name = phase_name - self.distance_trigger_units = distance_trigger_units self.event_channel_names = [ - # Dynamic.Mission.DISTANCE, Dynamic.Mission.MASS, ] self.num_events = len(self.event_channel_names) def event_equation_function(self, t, x): self.output_equation_function(t, x) - # 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() current_mass = self.get_val(Dynamic.Mission.MASS, units="lbm").squeeze() mass_trigger = 150000 @@ -94,7 +82,6 @@ def event_equation_function(self, t, x): alt_trigger = 50 return np.array([ current_alt - alt_trigger - # maybe mach ]) @@ -222,13 +209,8 @@ def test_phase(phases, ode_args_tab): for key, val in initial_values_cruise.items(): prob.set_val(key, **val) - # prob.set_val("traj.velocity", val=472, units="kn") - # prob.set_val("traj.velocity_rate", val=0, units="m/s**2") - # try: prob.run_model() - # except: - # prob.final_setup() with open('input_list.txt', 'w') as outfile: prob.model.list_inputs(out_stream=outfile,) @@ -267,16 +249,29 @@ def test_phase(phases, ode_args_tab): ode_args_tab['subsystem_options'] = {'core_aerodynamics': {'method': 'computed'}} # ode_args_tab['friction_key'] = Mission.Takeoff.ROLLING_FRICTION_COEFFICIENT # ode_args_tab['friction_key'] = Mission.Takeoff.BRAKING_FRICTION_COEFFICIENT + phase_kwargs = dict( ode_args=ode_args_tab, simupy_args=dict( DEBUG=True, ), ) + + def func(self, t, x): + self.output_equation_function(t, x) + current_mass = self.get_val(Dynamic.Mission.MASS, units="lbm").squeeze() + mass_trigger = 160000 + return np.array([ + current_mass - mass_trigger + ]) + + SGMCruise = SGMHeightEnergy + SGMCruise.event_equation_function = func + phase_vals = { } phases = {'HE': { - 'ode': SGMHeightEnergy(**phase_kwargs), + 'ode': SGMCruise(**phase_kwargs), 'vals_to_set': phase_vals }} diff --git a/aviary/mission/ode/specific_energy_rate.py b/aviary/mission/ode/specific_energy_rate.py index 62f0c5cf3..41002d0df 100644 --- a/aviary/mission/ode/specific_energy_rate.py +++ b/aviary/mission/ode/specific_energy_rate.py @@ -2,7 +2,6 @@ import openmdao.api as om from aviary.constants import GRAV_METRIC_FLOPS as gravity from aviary.variable_info.variables import Dynamic -from aviary.variable_info.enums import AnalysisScheme class SpecificEnergyRate(om.ExplicitComponent): @@ -13,8 +12,6 @@ class SpecificEnergyRate(om.ExplicitComponent): def initialize(self): self.options.declare('num_nodes', types=int) - self.options.declare("analysis_scheme", types=AnalysisScheme, default=AnalysisScheme.COLLOCATION, - desc="The analysis method that will be used to close the trajectory; for example collocation or time integration") def setup(self): nn = self.options['num_nodes'] From 62c1d20ec9b955fee618e3667a4c0a1c0cf9d2c2 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Mon, 5 Feb 2024 15:35:38 -0800 Subject: [PATCH 17/34] Adding flexible triggers and removing unused trajectories --- .../ode/time_integration_base_classes.py | 49 ++- .../phases/time_integration_traj.py | 289 ------------------ 2 files changed, 48 insertions(+), 290 deletions(-) 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 6b3a59afc..8179d6cac 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -27,10 +27,25 @@ def add_SGM_required_outputs(group: om.Group, outputs_to_add: dict): promotes_outputs=list(outputs_to_add.keys()), ) -# Subproblem used as a basis for forward in time integration phases. + +class event_trigger(): + def __init__( + self, + state: str, + value: 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, @@ -43,6 +58,7 @@ def __init__( parameters=None, outputs=None, controls=None, + triggers=None, include_state_outputs=False, rate_suffix="_rate", DEBUG=False, @@ -85,6 +101,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} @@ -343,6 +366,30 @@ 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'] + 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 event_equation_function(self, t, x): + self.output_equation_function(t, x) + event_values = [self.evaluate_trigger(trigger) for trigger in self.triggers] + return np.array(event_values) + + def evaluate_trigger(self, trigger: event_trigger): + trigger_value = trigger.value + if isinstance(trigger_value, str): + trigger_value = self.get_val(trigger_value, units=trigger.units).squeeze() + current_value = self.get_val(trigger.state, units=trigger.units) + return current_value - trigger_value + @property def get_val(self): return self.prob.get_val diff --git a/aviary/mission/gasp_based/phases/time_integration_traj.py b/aviary/mission/gasp_based/phases/time_integration_traj.py index 5e1884012..a58cb3367 100644 --- a/aviary/mission/gasp_based/phases/time_integration_traj.py +++ b/aviary/mission/gasp_based/phases/time_integration_traj.py @@ -104,292 +104,3 @@ def compute(self, inputs, outputs): 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 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 DEBUG: -# print("starting ascent") -# next_problem = ascent -# elif current_problem is ascent: -# next_problem = accel -# if DEBUG: -# print("starting accel") -# elif current_problem is accel: -# if DEBUG: -# print("climb1") -# next_problem = climb1 -# elif current_problem is climb1: -# if DEBUG: -# print("climb2") -# if event_idx != 0: -# if DEBUG: -# print("expected to hit alt trigger for climb1") -# next_problem = climb2 -# elif current_problem is climb2: -# if DEBUG: -# print("climb3") -# if event_idx != 1: -# if DEBUG: -# print("expected to hit speed trigger for climb2") -# next_problem = climb3 -# elif current_problem is climb3: -# if DEBUG: -# print("climb ending") -# if event_idx != 0: -# if DEBUG: -# print("expected to hit alt trigger for climb3") -# next_problem = None -# else: -# if 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 DEBUG: -# print("desc2") -# if event_idx != 1: -# if DEBUG: -# print("expected to hit speed trigger for desc1") -# sim_gen.send(desc2) -# elif current_problem is desc2: -# if DEBUG: -# print("desc3") -# if event_idx != 0: -# if DEBUG: -# print("expected to hit alt trigger for desc2") -# sim_gen.send(desc3) -# elif current_problem is desc3: -# if DEBUG: -# print("desc ending") -# if event_idx != 0: -# if DEBUG: -# print("expected to hit alt trigger for desc3") -# sim_gen.close() -# else: -# if DEBUG: -# print("unexpected termination") -# sim_gen.close() From a9e8b684f2620546612004af1975c360e5589669 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 6 Feb 2024 10:16:46 -0800 Subject: [PATCH 18/34] Using flexible triggers Included example of changing trigger value after creation --- .../phases/time_integration_phases.py | 54 +-------- .../ode/time_integration_base_classes.py | 9 +- .../phases/time_integration_phases.py | 105 ++++-------------- 3 files changed, 34 insertions(+), 134 deletions(-) diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index d7a9eb9f0..01683dd97 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -35,19 +35,7 @@ def __init__( **simupy_args) self.phase_name = phase_name - self.event_channel_names = [ - Dynamic.Mission.MASS, - ] - 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 = 150000 - return np.array([ - current_mass - mass_trigger - ]) + self.add_trigger(Dynamic.Mission.MASS, 150000, units='lbm') class SGMDetailedTakeoff(SimuPyProblem): @@ -71,18 +59,7 @@ def __init__( **simupy_args) self.phase_name = phase_name - self.event_channel_names = [ - Dynamic.Mission.ALTITUDE, - ] - self.num_events = len(self.event_channel_names) - - def event_equation_function(self, t, x): - self.output_equation_function(t, x) - current_alt = self.get_val(Dynamic.Mission.ALTITUDE, units="ft").squeeze() - alt_trigger = 50 - return np.array([ - current_alt - alt_trigger - ]) + self.add_trigger(Dynamic.Mission.ALTITUDE, 50, units='ft') class SGMDetailedLanding(SimuPyProblem): @@ -106,18 +83,7 @@ def __init__( **simupy_args) self.phase_name = phase_name - self.event_channel_names = [ - Dynamic.Mission.ALTITUDE, - ] - self.num_events = len(self.event_channel_names) - - def event_equation_function(self, t, x): - self.output_equation_function(t, x) - current_alt = self.get_val(Dynamic.Mission.ALTITUDE, units="ft").squeeze() - alt_trigger = 0 - return np.array([ - current_alt - alt_trigger - ]) + self.add_trigger(Dynamic.Mission.ALTITUDE, 0, units='ft') def test_phase(phases, ode_args_tab): @@ -257,21 +223,13 @@ def test_phase(phases, ode_args_tab): ), ) - def func(self, t, x): - self.output_equation_function(t, x) - current_mass = self.get_val(Dynamic.Mission.MASS, units="lbm").squeeze() - mass_trigger = 160000 - return np.array([ - current_mass - mass_trigger - ]) - - SGMCruise = SGMHeightEnergy - SGMCruise.event_equation_function = func + SGMCruise = SGMHeightEnergy(**phase_kwargs) + SGMCruise.triggers[0].value = 160000 phase_vals = { } phases = {'HE': { - 'ode': SGMCruise(**phase_kwargs), + 'ode': SGMCruise, 'vals_to_set': phase_vals }} 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 8179d6cac..3fd7cf032 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -378,6 +378,11 @@ def add_trigger(self, state, value, units=None, channel_name=None): 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 = [] + def event_equation_function(self, t, x): self.output_equation_function(t, x) event_values = [self.evaluate_trigger(trigger) for trigger in self.triggers] @@ -385,7 +390,9 @@ def event_equation_function(self, t, x): 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) + elif isinstance(trigger_value, str): trigger_value = self.get_val(trigger_value, units=trigger.units).squeeze() current_value = self.get_val(trigger.state, units=trigger.units) return current_value - trigger_value diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index a972070cd..5f46a4de8 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -50,15 +50,14 @@ def __init__( self.phase_name = phase_name self.VR_value = VR_value - self.VR_units = VR_units - self.event_channel_names = ["TAS"] - self.num_events = len(self.event_channel_names) + # self.VR_units = VR_units + self.add_trigger("TAS", "VR_value", units='ft/s') - def event_equation_function(self, t, x): - self.time = t - self.state = x - return self.get_val("TAS", units='ft/s') - self.VR_value - return self.get_val("TAS", units=self.VR_units) - self.VR_value + # def event_equation_function(self, t, x): + # self.time = t + # self.state = x + # return self.get_val("TAS", units='ft/s') - self.VR_value + # return self.get_val("TAS", units=self.VR_units) - self.VR_value class SGMRotation(SimuPyProblem): @@ -89,15 +88,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 @@ -389,16 +380,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): @@ -459,25 +441,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=self.speed_trigger_units) class SGMCruise(SimuPyProblem): @@ -528,26 +495,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): @@ -607,22 +557,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) From 28f9a786c9b6f38f8c10821c2b01de8aad6a3c9f Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 6 Feb 2024 10:47:24 -0800 Subject: [PATCH 19/34] fixed missing squeeze --- aviary/mission/gasp_based/ode/time_integration_base_classes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 3fd7cf032..04d1a76df 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -394,7 +394,7 @@ def evaluate_trigger(self, trigger: event_trigger): trigger_value = getattr(self, trigger_value) elif isinstance(trigger_value, str): trigger_value = self.get_val(trigger_value, units=trigger.units).squeeze() - current_value = self.get_val(trigger.state, units=trigger.units) + current_value = self.get_val(trigger.state, units=trigger.units).squeeze() return current_value - trigger_value @property From 52ea87ef87c2ab11f0b7216157a9c554301489df Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 6 Feb 2024 13:31:20 -0800 Subject: [PATCH 20/34] fixing units and type hint --- aviary/interface/default_phase_info/two_dof_fiti.py | 4 ++-- .../mission/gasp_based/ode/time_integration_base_classes.py | 4 +++- aviary/mission/gasp_based/phases/time_integration_phases.py | 2 +- 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/aviary/interface/default_phase_info/two_dof_fiti.py b/aviary/interface/default_phase_info/two_dof_fiti.py index 3c75a6d36..9e24d660b 100644 --- a/aviary/interface/default_phase_info/two_dof_fiti.py +++ b/aviary/interface/default_phase_info/two_dof_fiti.py @@ -65,7 +65,7 @@ def create_2dof_based_ascent_phases( climb1_vals = { 'alt_trigger': {'val': 10000, 'units': 'ft'}, 'EAS': {'val': 250, 'units': climb1_kwargs['input_speed_units']}, - 'speed_trigger': {'val': cruise_mach, 'units': None}, + 'speed_trigger': {'val': cruise_mach, 'units': 'unitless'}, } climb2_kwargs = dict( @@ -79,7 +79,7 @@ def create_2dof_based_ascent_phases( climb2_vals = { 'alt_trigger': {'val': cruise_alt, 'units': 'ft'}, 'EAS': {'val': 270, 'units': climb2_kwargs['input_speed_units']}, - 'speed_trigger': {'val': cruise_mach, 'units': None}, + 'speed_trigger': {'val': cruise_mach, 'units': 'unitless'}, } climb3_kwargs = dict( 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 04d1a76df..a016bcd81 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -32,7 +32,7 @@ class event_trigger(): def __init__( self, state: str, - value: str | float | int, + value: tuple[str, float, int], units: str, channel_name: str = None, ): @@ -369,6 +369,8 @@ def update_equation_function(self, t, x, event_channels=None): 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 diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index efeef26c9..f2637c771 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -443,7 +443,7 @@ def __init__( 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) + units="speed_trigger_units") class SGMCruise(SimuPyProblem): From 6042876307d10d3c5070e8f712199a1e7a0ae1d3 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Thu, 8 Feb 2024 10:06:03 -0800 Subject: [PATCH 21/34] Added some temporary debugging statements and tweaked triggers --- .../default_phase_info/two_dof_fiti.py | 5 ++-- aviary/interface/methods_for_level2.py | 20 ++++++++------- .../phases/time_integration_phases.py | 3 +++ aviary/mission/gasp_based/ode/ascent_eom.py | 16 ++++++++---- aviary/mission/gasp_based/ode/ascent_ode.py | 18 +++++++++++++ .../mission/gasp_based/ode/groundroll_ode.py | 20 +++++++++++++++ .../ode/time_integration_base_classes.py | 24 ++++++++++++------ .../phases/time_integration_phases.py | 25 +++++++++++++++++-- .../phases/time_integration_traj.py | 23 +++++++++++++++-- 9 files changed, 126 insertions(+), 28 deletions(-) diff --git a/aviary/interface/default_phase_info/two_dof_fiti.py b/aviary/interface/default_phase_info/two_dof_fiti.py index 9e24d660b..9f312003c 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', @@ -65,7 +66,7 @@ def create_2dof_based_ascent_phases( climb1_vals = { 'alt_trigger': {'val': 10000, 'units': 'ft'}, 'EAS': {'val': 250, 'units': climb1_kwargs['input_speed_units']}, - 'speed_trigger': {'val': cruise_mach, 'units': 'unitless'}, + 'speed_trigger': {'val': cruise_mach, 'units': None}, } climb2_kwargs = dict( @@ -79,7 +80,7 @@ def create_2dof_based_ascent_phases( climb2_vals = { 'alt_trigger': {'val': cruise_alt, 'units': 'ft'}, 'EAS': {'val': 270, 'units': climb2_kwargs['input_speed_units']}, - 'speed_trigger': {'val': cruise_mach, 'units': 'unitless'}, + 'speed_trigger': {'val': cruise_mach, 'units': None}, } climb3_kwargs = dict( diff --git a/aviary/interface/methods_for_level2.py b/aviary/interface/methods_for_level2.py index 8d2c32eb7..2c5c55f5f 100644 --- a/aviary/interface/methods_for_level2.py +++ b/aviary/interface/methods_for_level2.py @@ -216,7 +216,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 @@ -929,16 +929,18 @@ def add_phases(self, phase_info_parameterization=None): self.ode_args, cruise_mach=self.cruise_mach) - descent_estimation = descent_range_and_fuel( - phases=descent_phases, - initial_mass=initial_mass, - cruise_alt=self.cruise_alt, - cruise_mach=self.cruise_mach) + # TEMPORARY + # descent_estimation = descent_range_and_fuel( + # phases=descent_phases, + # initial_mass=initial_mass, + # cruise_alt=self.cruise_alt, + # cruise_mach=self.cruise_mach) - estimated_descent_range = descent_estimation['refined_guess']['distance_flown'] - end_of_cruise_range = self.target_range - estimated_descent_range + # estimated_descent_range = descent_estimation['refined_guess']['distance_flown'] + # end_of_cruise_range = self.target_range - estimated_descent_range - estimated_descent_fuel = descent_estimation['refined_guess']['fuel_burned'] + # estimated_descent_fuel = descent_estimation['refined_guess']['fuel_burned'] + estimated_descent_fuel = 250 cruise_kwargs = dict( input_speed_type=SpeedType.MACH, diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index 01683dd97..6f8edafaf 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -24,6 +24,7 @@ def __init__( super().__init__(MissionODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), + problem_name=phase_name, outputs=[], states=[ Dynamic.Mission.MASS, @@ -48,6 +49,7 @@ def __init__( super().__init__(TakeoffODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), + problem_name=phase_name, outputs=[], states=[ Dynamic.Mission.MASS, @@ -72,6 +74,7 @@ def __init__( super().__init__(LandingODE( analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), + problem_name=phase_name, outputs=[], states=[ Dynamic.Mission.MASS, diff --git a/aviary/mission/gasp_based/ode/ascent_eom.py b/aviary/mission/gasp_based/ode/ascent_eom.py index ea6f7c3cd..3f89b8ade 100644 --- a/aviary/mission/gasp_based/ode/ascent_eom.py +++ b/aviary/mission/gasp_based/ode/ascent_eom.py @@ -151,11 +151,17 @@ 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 - / (TAS * weight) - ) + import warnings + warnings.simplefilter("error") + try: + outputs[Dynamic.Mission.FLIGHT_PATH_ANGLE_RATE] = ( + (thrust_across_flightpath + incremented_lift - weight * np.cos(gamma)) + * GRAV_ENGLISH_GASP + / (TAS * weight) + ) + except: + print(TAS, weight) + exit() outputs[Dynamic.Mission.ALTITUDE_RATE] = TAS * np.sin(gamma) outputs[Dynamic.Mission.DISTANCE_RATE] = TAS * np.cos(gamma) diff --git a/aviary/mission/gasp_based/ode/ascent_ode.py b/aviary/mission/gasp_based/ode/ascent_ode.py index 14cfab863..2d4a4e062 100644 --- a/aviary/mission/gasp_based/ode/ascent_ode.py +++ b/aviary/mission/gasp_based/ode/ascent_ode.py @@ -108,6 +108,24 @@ def setup(self): self.add_excess_rate_comps(nn) + if analysis_scheme is AnalysisScheme.SHOOTING or False: + from aviary.utils.functions import create_printcomp + dummy_comp = create_printcomp( + all_inputs=[ + "t_curr", + 'alpha', + Dynamic.Mission.VELOCITY, + Dynamic.Mission.ALTITUDE, + ], + input_units={ + 't_curr': 's', + 'alpha': 'deg', + }) + self.add_subsystem( + "dummy_comp", + dummy_comp(), + promotes_inputs=["*"],) + ParamPort.set_default_vals(self) self.set_input_defaults("t_init_flaps", val=47.5) self.set_input_defaults("t_init_gear", val=37.3) diff --git a/aviary/mission/gasp_based/ode/groundroll_ode.py b/aviary/mission/gasp_based/ode/groundroll_ode.py index bb9700dea..f9c8d0e3b 100644 --- a/aviary/mission/gasp_based/ode/groundroll_ode.py +++ b/aviary/mission/gasp_based/ode/groundroll_ode.py @@ -103,6 +103,26 @@ def setup(self): Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL), "dt_dv"]) + import aviary.api as av + if analysis_scheme is av.AnalysisScheme.SHOOTING and False: + from aviary.utils.functions import create_printcomp + dummy_comp = create_printcomp( + all_inputs=[ + "t_curr", + 'alpha', + 'velocity' + ], + input_units={ + 't_curr': 's', + 'alpha': 'deg', + }) + self.add_subsystem( + "dummy_comp", + dummy_comp(), + promotes_inputs=["*"],) + # self.set_input_defaults('start_rotation', val=0, units='s') + # self.set_input_defaults('rotation_rate', val=10/3, units='deg/s') + ParamPort.set_default_vals(self) self.set_input_defaults("t_init_flaps", val=100.) self.set_input_defaults("t_init_gear", val=100.) 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 a016bcd81..c8b8d408b 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -49,6 +49,7 @@ class SimuPyProblem(SimulationMixin): def __init__( self, ode, + problem_name='', time_independent=False, t_name="t_curr", states=None, @@ -220,11 +221,14 @@ def __init__( # TODO: add defensive checks to make sure dimensions match in both setup and # calls - if DEBUG: - om.n2(prob, outfile="n2_simupy_problem.html", show_browser=False) - with open('input_list_simupy.txt', 'w') as outfile: + if DEBUG or True: + 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) + # print(states) @property def time(self): @@ -387,15 +391,19 @@ def clear_triggers(self): 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 hasattr(self, trigger_value): - trigger_value = getattr(self, trigger_value) - elif isinstance(trigger_value, str): - trigger_value = self.get_val(trigger_value, units=trigger.units).squeeze() + 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 diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index f2637c771..743f9bbcf 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -35,6 +35,7 @@ def __init__( super().__init__( GroundrollODE(analysis_scheme=AnalysisScheme.SHOOTING, **ode_args), + problem_name=phase_name, outputs=["normal_force"], states=[ Dynamic.Mission.MASS, @@ -51,7 +52,9 @@ def __init__( self.phase_name = phase_name self.VR_value = VR_value # self.VR_units = VR_units - self.add_trigger(Dynamic.Mission.VELOCITY, "VR_value", units='ft/s') + 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 @@ -74,11 +77,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={ @@ -88,6 +93,15 @@ def __init__( self.phase_name = phase_name self.add_trigger("normal_force", 0, units='lbf') + # 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 # TODO : turn these into parameters? inputs? they need to match between @@ -119,6 +133,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", @@ -129,6 +144,7 @@ def __init__( Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, + Dynamic.Mission.VELOCITY, ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ @@ -195,7 +211,8 @@ 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 @@ -367,6 +384,7 @@ 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, @@ -417,6 +435,7 @@ def __init__( self.alt_trigger_units = alt_trigger_units super().__init__( ode, + problem_name=phase_name, outputs=[ "alpha", Dynamic.Mission.FLIGHT_PATH_ANGLE, @@ -473,6 +492,7 @@ def __init__( super().__init__( ode, + problem_name=phase_name, outputs=[ "alpha", # ? "lift", @@ -534,6 +554,7 @@ def __init__( self.alt_trigger_units = alt_trigger_units super().__init__( ode, + problem_name=phase_name, outputs=[ "alpha", "required_lift", diff --git a/aviary/mission/gasp_based/phases/time_integration_traj.py b/aviary/mission/gasp_based/phases/time_integration_traj.py index a58cb3367..e2238dc7e 100644 --- a/aviary/mission/gasp_based/phases/time_integration_traj.py +++ b/aviary/mission/gasp_based/phases/time_integration_traj.py @@ -7,7 +7,7 @@ from aviary.utils.aviary_values import AviaryValues from aviary.variable_info.variables import Aircraft, Mission, Dynamic -DEBUG = 0 +DEBUG = 1 class TimeIntegrationTrajBase(SGMTrajBase): @@ -70,7 +70,12 @@ def compute(self, inputs, outputs): phase.rotation.set_val(name.replace( 'rotation.', ''), data['val'], units=data['units']) else: - phase.set_val(name, data['val'], units=data['units']) + try: + phase.set_val(name, data['val'], units=data['units']) + except TypeError: + print('\n') + print(phase_name, name, data['val'], data['units']) + exit() ode_index = 0 sim_gen = self.compute_traj_loop(self.ODEs[0], inputs, outputs) @@ -89,12 +94,26 @@ def compute(self, inputs, outputs): next_problem = None print('Finished: '+current_problem.phase_name) + print([name+'_'+data['units'] + for name, data in current_problem.states.items()]) + print(x_final) + print('outputs') + for output_name, unit in current_problem.outputs.items(): + val = current_problem.get_val(output_name, units=unit)[0] + print(output_name+':', val, unit) + + # print(self.outputs.items()) + # print(current_problem.output) if next_problem is not None: if type(current_problem) is SGMGroundroll: next_problem.prob.set_val("start_rotation", t_start_rotation) + # next_problem.prob.set_val(Dynamic.Mission.VELOCITY, 150, units='kn') elif type(current_problem) is SGMRotation: next_problem.rotation.set_val("start_rotation", t_start_rotation) + print('\n\n') print('Starting: '+next_problem.phase_name) + print([name+'_'+data['units'] + for name, data in next_problem.states.items()]) sim_gen.send(next_problem) else: print('Reached the end of the Trajectory!') From 83352702322b4b14a9ac5c4cd5f79563168e3636 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 9 Feb 2024 17:12:15 -0800 Subject: [PATCH 22/34] added debug environment files to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) 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 From 33c0e8bec3b24a66ad5f0f194c14decc5cbe04d9 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Mon, 12 Feb 2024 14:55:43 -0800 Subject: [PATCH 23/34] fixed typo --- aviary/mission/gasp_based/ode/time_integration_base_classes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 c8b8d408b..62b2cd1b0 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -161,7 +161,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(): From 35bc9dd492df6099176947f4a0740fd921d4b5ad Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 13 Feb 2024 16:29:44 -0800 Subject: [PATCH 24/34] adding phase names to ascent subphases and adding states to phases --- aviary/mission/gasp_based/ode/ascent_eom.py | 17 ++++++----------- aviary/mission/gasp_based/ode/ascent_ode.py | 18 ------------------ aviary/mission/gasp_based/ode/climb_ode.py | 2 -- aviary/mission/gasp_based/ode/descent_ode.py | 2 -- .../mission/gasp_based/ode/flight_path_ode.py | 3 --- .../phases/time_integration_phases.py | 16 ++++++++++------ 6 files changed, 16 insertions(+), 42 deletions(-) diff --git a/aviary/mission/gasp_based/ode/ascent_eom.py b/aviary/mission/gasp_based/ode/ascent_eom.py index 3f89b8ade..0cda876c5 100644 --- a/aviary/mission/gasp_based/ode/ascent_eom.py +++ b/aviary/mission/gasp_based/ode/ascent_eom.py @@ -151,17 +151,12 @@ def compute(self, inputs, outputs): * GRAV_ENGLISH_GASP / weight ) - import warnings - warnings.simplefilter("error") - try: - outputs[Dynamic.Mission.FLIGHT_PATH_ANGLE_RATE] = ( - (thrust_across_flightpath + incremented_lift - weight * np.cos(gamma)) - * GRAV_ENGLISH_GASP - / (TAS * weight) - ) - except: - print(TAS, weight) - exit() + + outputs[Dynamic.Mission.FLIGHT_PATH_ANGLE_RATE] = ( + (thrust_across_flightpath + incremented_lift - weight * np.cos(gamma)) + * GRAV_ENGLISH_GASP + / (TAS * weight) + ) outputs[Dynamic.Mission.ALTITUDE_RATE] = TAS * np.sin(gamma) outputs[Dynamic.Mission.DISTANCE_RATE] = TAS * np.cos(gamma) diff --git a/aviary/mission/gasp_based/ode/ascent_ode.py b/aviary/mission/gasp_based/ode/ascent_ode.py index 2d4a4e062..14cfab863 100644 --- a/aviary/mission/gasp_based/ode/ascent_ode.py +++ b/aviary/mission/gasp_based/ode/ascent_ode.py @@ -108,24 +108,6 @@ def setup(self): self.add_excess_rate_comps(nn) - if analysis_scheme is AnalysisScheme.SHOOTING or False: - from aviary.utils.functions import create_printcomp - dummy_comp = create_printcomp( - all_inputs=[ - "t_curr", - 'alpha', - Dynamic.Mission.VELOCITY, - Dynamic.Mission.ALTITUDE, - ], - input_units={ - 't_curr': 's', - 'alpha': 'deg', - }) - self.add_subsystem( - "dummy_comp", - dummy_comp(), - promotes_inputs=["*"],) - ParamPort.set_default_vals(self) self.set_input_defaults("t_init_flaps", val=47.5) self.set_input_defaults("t_init_gear", val=37.3) 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/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index 743f9bbcf..f06ee7b0a 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -145,6 +145,8 @@ def __init__( Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, Dynamic.Mission.VELOCITY, + Dynamic.Mission.FLIGHT_PATH_ANGLE, + "alpha", ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ @@ -217,15 +219,15 @@ def __init__( 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,) @@ -390,6 +392,7 @@ def __init__( Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, + Dynamic.Mission.VELOCITY, ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ @@ -506,6 +509,7 @@ def __init__( Dynamic.Mission.MASS, Dynamic.Mission.DISTANCE, Dynamic.Mission.ALTITUDE, + Dynamic.Mission.VELOCITY, ], # state_units=['lbm','nmi','ft'], alternate_state_rate_names={ From edec53741c7870cc70cbdfead6fb5c95d00530ef Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 13 Feb 2024 16:30:46 -0800 Subject: [PATCH 25/34] restoring descent estimation --- aviary/interface/methods_for_level2.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/aviary/interface/methods_for_level2.py b/aviary/interface/methods_for_level2.py index 2c5c55f5f..2bf28531d 100644 --- a/aviary/interface/methods_for_level2.py +++ b/aviary/interface/methods_for_level2.py @@ -929,17 +929,16 @@ def add_phases(self, phase_info_parameterization=None): self.ode_args, cruise_mach=self.cruise_mach) - # TEMPORARY - # descent_estimation = descent_range_and_fuel( - # phases=descent_phases, - # initial_mass=initial_mass, - # cruise_alt=self.cruise_alt, - # cruise_mach=self.cruise_mach) + descent_estimation = descent_range_and_fuel( + phases=descent_phases, + initial_mass=initial_mass, + cruise_alt=self.cruise_alt, + cruise_mach=self.cruise_mach) - # estimated_descent_range = descent_estimation['refined_guess']['distance_flown'] - # end_of_cruise_range = self.target_range - estimated_descent_range + estimated_descent_range = descent_estimation['refined_guess']['distance_flown'] + end_of_cruise_range = self.target_range - estimated_descent_range - # estimated_descent_fuel = descent_estimation['refined_guess']['fuel_burned'] + estimated_descent_fuel = descent_estimation['refined_guess']['fuel_burned'] estimated_descent_fuel = 250 cruise_kwargs = dict( From 546f1fa2864ff002c739ce05879065020d942f6e Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 13 Feb 2024 16:42:14 -0800 Subject: [PATCH 26/34] reverted temporary debugging --- aviary/interface/methods_for_level2.py | 1 - .../mission/gasp_based/ode/groundroll_ode.py | 20 ------------- .../ode/time_integration_base_classes.py | 4 +-- .../phases/time_integration_phases.py | 9 ------ .../phases/time_integration_traj.py | 30 ++++++++----------- 5 files changed, 14 insertions(+), 50 deletions(-) diff --git a/aviary/interface/methods_for_level2.py b/aviary/interface/methods_for_level2.py index 2bf28531d..f5140eea5 100644 --- a/aviary/interface/methods_for_level2.py +++ b/aviary/interface/methods_for_level2.py @@ -939,7 +939,6 @@ def add_phases(self, phase_info_parameterization=None): end_of_cruise_range = self.target_range - estimated_descent_range estimated_descent_fuel = descent_estimation['refined_guess']['fuel_burned'] - estimated_descent_fuel = 250 cruise_kwargs = dict( input_speed_type=SpeedType.MACH, diff --git a/aviary/mission/gasp_based/ode/groundroll_ode.py b/aviary/mission/gasp_based/ode/groundroll_ode.py index f9c8d0e3b..bb9700dea 100644 --- a/aviary/mission/gasp_based/ode/groundroll_ode.py +++ b/aviary/mission/gasp_based/ode/groundroll_ode.py @@ -103,26 +103,6 @@ def setup(self): Dynamic.Mission.FUEL_FLOW_RATE_NEGATIVE_TOTAL), "dt_dv"]) - import aviary.api as av - if analysis_scheme is av.AnalysisScheme.SHOOTING and False: - from aviary.utils.functions import create_printcomp - dummy_comp = create_printcomp( - all_inputs=[ - "t_curr", - 'alpha', - 'velocity' - ], - input_units={ - 't_curr': 's', - 'alpha': 'deg', - }) - self.add_subsystem( - "dummy_comp", - dummy_comp(), - promotes_inputs=["*"],) - # self.set_input_defaults('start_rotation', val=0, units='s') - # self.set_input_defaults('rotation_rate', val=10/3, units='deg/s') - ParamPort.set_default_vals(self) self.set_input_defaults("t_init_flaps", val=100.) self.set_input_defaults("t_init_gear", val=100.) 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 62b2cd1b0..60b03f598 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -221,14 +221,14 @@ def __init__( # TODO: add defensive checks to make sure dimensions match in both setup and # calls - if DEBUG or True: + if DEBUG: 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) + print(states) @property def time(self): diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index f06ee7b0a..234662ded 100644 --- a/aviary/mission/gasp_based/phases/time_integration_phases.py +++ b/aviary/mission/gasp_based/phases/time_integration_phases.py @@ -93,15 +93,6 @@ def __init__( self.phase_name = phase_name self.add_trigger("normal_force", 0, units='lbf') - # 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 # TODO : turn these into parameters? inputs? they need to match between diff --git a/aviary/mission/gasp_based/phases/time_integration_traj.py b/aviary/mission/gasp_based/phases/time_integration_traj.py index e2238dc7e..07fea5cf1 100644 --- a/aviary/mission/gasp_based/phases/time_integration_traj.py +++ b/aviary/mission/gasp_based/phases/time_integration_traj.py @@ -7,7 +7,7 @@ from aviary.utils.aviary_values import AviaryValues from aviary.variable_info.variables import Aircraft, Mission, Dynamic -DEBUG = 1 +DEBUG = 0 class TimeIntegrationTrajBase(SGMTrajBase): @@ -70,12 +70,7 @@ def compute(self, inputs, outputs): phase.rotation.set_val(name.replace( 'rotation.', ''), data['val'], units=data['units']) else: - try: - phase.set_val(name, data['val'], units=data['units']) - except TypeError: - print('\n') - print(phase_name, name, data['val'], data['units']) - exit() + phase.set_val(name, data['val'], units=data['units']) ode_index = 0 sim_gen = self.compute_traj_loop(self.ODEs[0], inputs, outputs) @@ -94,26 +89,25 @@ def compute(self, inputs, outputs): next_problem = None print('Finished: '+current_problem.phase_name) - print([name+'_'+data['units'] - for name, data in current_problem.states.items()]) - print(x_final) - print('outputs') - for output_name, unit in current_problem.outputs.items(): - val = current_problem.get_val(output_name, units=unit)[0] - print(output_name+':', val, unit) + # print([name+'_'+data['units'] + # for name, data in current_problem.states.items()]) + # print(x_final) + # print('outputs') + # for output_name, unit in current_problem.outputs.items(): + # val = current_problem.get_val(output_name, units=unit)[0] + # print(output_name+':', val, unit) # print(self.outputs.items()) # print(current_problem.output) if next_problem is not None: if type(current_problem) is SGMGroundroll: next_problem.prob.set_val("start_rotation", t_start_rotation) - # next_problem.prob.set_val(Dynamic.Mission.VELOCITY, 150, units='kn') elif type(current_problem) is SGMRotation: next_problem.rotation.set_val("start_rotation", t_start_rotation) - print('\n\n') + # print('\n\n') print('Starting: '+next_problem.phase_name) - print([name+'_'+data['units'] - for name, data in next_problem.states.items()]) + # print([name+'_'+data['units'] + # for name, data in next_problem.states.items()]) sim_gen.send(next_problem) else: print('Reached the end of the Trajectory!') From 3ec4f9bb82016e33d3cfc63045884592357d8fd4 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Wed, 14 Feb 2024 09:40:00 -0800 Subject: [PATCH 27/34] updated debug to verbosity --- .../mission/flops_based/phases/time_integration_phases.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/aviary/mission/flops_based/phases/time_integration_phases.py b/aviary/mission/flops_based/phases/time_integration_phases.py index 6f8edafaf..5c02e89c9 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -6,7 +6,7 @@ from aviary.mission.flops_based.ode.takeoff_ode import TakeoffODE from aviary.mission.gasp_based.phases.time_integration_traj import FlexibleTraj from aviary.mission.gasp_based.ode.time_integration_base_classes import SimuPyProblem -from aviary.variable_info.enums import AlphaModes, AnalysisScheme, SpeedType +from aviary.variable_info.enums import AnalysisScheme, Verbosity from aviary.variable_info.variables import Aircraft, Dynamic, Mission from aviary.variable_info.variables_in import VariablesIn from aviary.subsystems.premission import CorePreMission @@ -203,7 +203,7 @@ def test_phase(phases, ode_args_tab): aviary_inputs, initial_guesses = create_vehicle( 'models/test_aircraft/aircraft_for_bench_FwFm.csv') - aviary_inputs.set_val('debug_mode', False) + # aviary_inputs.set_val('debug_mode', False) 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, @@ -222,7 +222,7 @@ def test_phase(phases, ode_args_tab): phase_kwargs = dict( ode_args=ode_args_tab, simupy_args=dict( - DEBUG=True, + verbosity=Verbosity.QUIET, ), ) From 085694035d423f3be2ac0b33a6a0682d5f3c155d Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Wed, 14 Feb 2024 14:31:15 -0800 Subject: [PATCH 28/34] added a bench test for SGM and add some key to list conversion --- .../ode/time_integration_base_classes.py | 8 +++---- .../benchmark_tests/test_bench_GwGm.py | 23 +++++++++++++++++++ 2 files changed, 27 insertions(+), 4 deletions(-) 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 dbd50dac4..67c0082ca 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -697,7 +697,7 @@ def compute_partials(self, inputs, J): next_prob.state_equation_function(next_res.t[-1], next_res.x[-1, :]) costate[:] = next_prob.compute_totals( output, - next_prob.states.keys(), + list(next_prob.states.keys()), return_format='array' ).squeeze() @@ -743,7 +743,7 @@ def compute_partials(self, inputs, J): else: dg_dx[0, :] = prob.compute_totals( [channel_name], - prob.states.keys(), + list(prob.states.keys()), return_format='array' ) @@ -804,7 +804,7 @@ def compute_partials(self, inputs, J): dh_j_dx = prob.compute_totals( [state_name], - prob.states.keys(), + list(prob.states.keys()), return_format='array').squeeze() dh_dparam[state_idx, :] = prob.compute_totals( @@ -835,7 +835,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(), + list(prob.states.keys()), return_format='array').T if param_dict: df_dparam_data[idx, ...] = prob.compute_totals( diff --git a/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py b/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py index db7bde97d..6ccf14ecf 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="SNOPT") + 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() From cfdb0daf3650249efc56d2d1e6c5585d5b120f16 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 16 Feb 2024 15:56:46 -0800 Subject: [PATCH 29/34] moving HE SGM example to test case and minor refactoring --- aviary/mission/flops_based/ode/mission_ODE.py | 30 --- .../test/test_time_integration_phases.py | 204 ++++++++++++++++++ .../phases/time_integration_phases.py | 173 +-------------- .../ode/time_integration_base_classes.py | 11 +- .../phases/time_integration_phases.py | 5 +- .../phases/time_integration_traj.py | 44 ++-- .../test/test_idle_descent_estimation.py | 4 +- 7 files changed, 242 insertions(+), 229 deletions(-) create mode 100644 aviary/mission/flops_based/phases/test/test_time_integration_phases.py diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index 6b8068fd6..7b97787ea 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -185,36 +185,6 @@ def setup(self): promotes_inputs=['*'], promotes_outputs=['*']) - self.set_input_defaults(Mission.Design.GROSS_MASS, val=1, units='kg') - self.set_input_defaults( - Aircraft.Fuselage.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.Fuselage.FINENESS, val=1, units='unitless') - self.set_input_defaults(Aircraft.Fuselage.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults( - Aircraft.VerticalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.VerticalTail.FINENESS, val=1, units='unitless') - self.set_input_defaults(Aircraft.VerticalTail.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults( - Aircraft.HorizontalTail.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.HorizontalTail.FINENESS, - val=1, units='unitless') - self.set_input_defaults( - Aircraft.HorizontalTail.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults(Aircraft.Wing.CHARACTERISTIC_LENGTH, val=1, units='ft') - self.set_input_defaults(Aircraft.Wing.FINENESS, val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.WETTED_AREA, val=1, units='ft**2') - self.set_input_defaults( - Aircraft.Wing.SPAN_EFFICIENCY_FACTOR, val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.TAPER_RATIO, val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.THICKNESS_TO_CHORD, - val=1, units='unitless') - self.set_input_defaults(Aircraft.Wing.SWEEP, val=1, units='deg') - self.set_input_defaults(Aircraft.Wing.ASPECT_RATIO, val=1, units='unitless') - self.set_input_defaults( - Aircraft.Design.LIFT_DEPENDENT_DRAG_COEFF_FACTOR, val=1, units='unitless') - self.set_input_defaults( - Aircraft.Design.ZERO_LIFT_DRAG_COEFF_FACTOR, val=1, units='unitless') - 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') 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..68af92190 --- /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 phase_info, 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 + + +class IdleDescentTestCase(unittest.TestCase): + def setUp(self): + core_subsystems = [prop, geom, aero] + + 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=core_subsystems) + engine = EngineDeck(options=aviary_inputs) + preprocess_propulsion(aviary_inputs, [engine]) + + 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': 0, 'units': "m/s"}, + # } + + # ode_args = self.ode_args + # ode_args['friction_key'] = Mission.Takeoff.ROLLING_FRICTION_COEFFICIENT + # increased_alt_Takeoff = SGMDetailedTakeoff( + # ode_args, + # simupy_args=dict(verbosity=Verbosity.QUIET,) + # ) + # increased_alt_Takeoff.clear_triggers() + # increased_alt_Takeoff.add_trigger(Dynamic.Mission.ALTITUDE, 500, units='ft') + + # phases = {'HE': { + # 'ode': increased_alt_Takeoff, + # 'vals_to_set': {} + # }} + + # final_states = self.run_simulation(phases, initial_values_takeoff) + # assert_near_equal(final_states['altitude'], 500, 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 index 5c02e89c9..7590d469c 100644 --- a/aviary/mission/flops_based/phases/time_integration_phases.py +++ b/aviary/mission/flops_based/phases/time_integration_phases.py @@ -1,24 +1,16 @@ -import numpy as np -import openmdao.api as om - 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.phases.time_integration_traj import FlexibleTraj from aviary.mission.gasp_based.ode.time_integration_base_classes import SimuPyProblem -from aviary.variable_info.enums import AnalysisScheme, Verbosity -from aviary.variable_info.variables import Aircraft, Dynamic, Mission -from aviary.variable_info.variables_in import VariablesIn -from aviary.subsystems.premission import CorePreMission -from aviary.utils.functions import set_aviary_initial_values -import warnings +from aviary.variable_info.enums import AnalysisScheme +from aviary.variable_info.variables import Dynamic class SGMHeightEnergy(SimuPyProblem): def __init__( self, - phase_name='cruise', - ode_args={}, + ode_args, + phase_name='mission', simupy_args={}, ): super().__init__(MissionODE( @@ -33,6 +25,7 @@ def __init__( ], 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 @@ -42,8 +35,8 @@ def __init__( class SGMDetailedTakeoff(SimuPyProblem): def __init__( self, + ode_args, phase_name='detailed_takeoff', - ode_args={}, simupy_args={}, ): super().__init__(TakeoffODE( @@ -58,6 +51,7 @@ def __init__( ], 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 @@ -67,8 +61,8 @@ def __init__( class SGMDetailedLanding(SimuPyProblem): def __init__( self, + ode_args, phase_name='detailed_landing', - ode_args={}, simupy_args={}, ): super().__init__(LandingODE( @@ -83,157 +77,8 @@ def __init__( ], 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') - - -def test_phase(phases, ode_args_tab): - 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 = ode_args_tab['aviary_options'] - subsystems = ode_args_tab['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 = om.Group() - 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=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, aviary_inputs, meta_data=BaseMetaData) - - warnings.simplefilter("ignore", om.PromotionWarning) - - prob.setup() - - 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': 0, 'units': "m/s"}, - } - 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"}, - } - 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"}, - } - - for key, val in initial_values_cruise.items(): - prob.set_val(key, **val) - - prob.run_model() - - with open('input_list.txt', 'w') as outfile: - prob.model.list_inputs(out_stream=outfile,) - with open('output_list.txt', 'w') as outfile: - prob.model.list_outputs(out_stream=outfile) - - final_distance = prob.get_val('traj.distance_final', units='NM')[0] - final_mass = prob.get_val('traj.mass_final', units='lbm')[0] - print(final_distance, final_mass) - - -if __name__ == '__main__': - from aviary.interface.default_phase_info.height_energy import phase_info, aero, prop, geom - from aviary.subsystems.propulsion.engine_deck import EngineDeck - from aviary.variable_info.variables import Aircraft, Mission, Dynamic - 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 - - core_subsystems = [prop, geom, aero] - - aviary_inputs, initial_guesses = create_vehicle( - 'models/test_aircraft/aircraft_for_bench_FwFm.csv') - # aviary_inputs.set_val('debug_mode', False) - 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_tab = dict(aviary_options=aviary_inputs, core_subsystems=core_subsystems) - engine = EngineDeck(options=aviary_inputs) - preprocess_propulsion(aviary_inputs, [engine]) - - ode_args_tab['num_nodes'] = 1 - ode_args_tab['subsystem_options'] = {'core_aerodynamics': {'method': 'computed'}} - # ode_args_tab['friction_key'] = Mission.Takeoff.ROLLING_FRICTION_COEFFICIENT - # ode_args_tab['friction_key'] = Mission.Takeoff.BRAKING_FRICTION_COEFFICIENT - - phase_kwargs = dict( - ode_args=ode_args_tab, - simupy_args=dict( - verbosity=Verbosity.QUIET, - ), - ) - - SGMCruise = SGMHeightEnergy(**phase_kwargs) - SGMCruise.triggers[0].value = 160000 - - phase_vals = { - } - phases = {'HE': { - 'ode': SGMCruise, - 'vals_to_set': phase_vals - }} - - test_phase(phases, ode_args_tab) 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 67c0082ca..47ae1ec5b 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -7,6 +7,7 @@ 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): @@ -50,6 +51,8 @@ class SimuPyProblem(SimulationMixin): def __init__( self, ode, + aviary_options=None, + meta_data=_MetaData, problem_name='', time_independent=False, t_name="t_curr", @@ -92,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, @@ -418,12 +425,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 diff --git a/aviary/mission/gasp_based/phases/time_integration_phases.py b/aviary/mission/gasp_based/phases/time_integration_phases.py index e0d42e66c..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 diff --git a/aviary/mission/gasp_based/phases/time_integration_traj.py b/aviary/mission/gasp_based/phases/time_integration_traj.py index 51d425d4c..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,32 +83,25 @@ def compute(self, inputs, outputs): except IndexError: next_problem = None - print('Finished: '+current_problem.phase_name) - # print([name+'_'+data['units'] - # for name, data in current_problem.states.items()]) - # print(x_final) - # print('outputs') - # for output_name, unit in current_problem.outputs.items(): - # val = current_problem.get_val(output_name, units=unit)[0] - # print(output_name+':', val, unit) - - # print(self.outputs.items()) - # print(current_problem.output) + 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('\n\n') - print('Starting: '+next_problem.phase_name) - # print([name+'_'+data['units'] - # for name, data in next_problem.states.items()]) + + 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) + if self.verbosity.value >= 1: + print('t_final', t_final) + print('x_final', x_final) + print(self.ODEs[-1].states) diff --git a/aviary/mission/gasp_based/test/test_idle_descent_estimation.py b/aviary/mission/gasp_based/test/test_idle_descent_estimation.py index c748d6d4b..d1363e8c9 100644 --- a/aviary/mission/gasp_based/test/test_idle_descent_estimation.py +++ b/aviary/mission/gasp_based/test/test_idle_descent_estimation.py @@ -29,8 +29,8 @@ def test_case1(self): results = descent_range_and_fuel(ode_args=ode_args)['refined_guess'] # Values obtained by running idle_descent_estimation - assert_near_equal(results['distance_flown'], 91.8911599691433, tol) - assert_near_equal(results['fuel_burned'], 236.73893823639082, tol) + assert_near_equal(results['distance_flown'], 91.50393932590772, tol) + assert_near_equal(results['fuel_burned'], 235.39517919329228, tol) if __name__ == "__main__": From 5e912d0e18895dc7f3021873c18ba097a389bbf0 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Fri, 16 Feb 2024 16:05:14 -0800 Subject: [PATCH 30/34] adding skip for missing pyoptsparse --- .../flops_based/phases/test/test_time_integration_phases.py | 2 ++ 1 file changed, 2 insertions(+) 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 index 68af92190..376db30fa 100644 --- a/aviary/mission/flops_based/phases/test/test_time_integration_phases.py +++ b/aviary/mission/flops_based/phases/test/test_time_integration_phases.py @@ -19,8 +19,10 @@ import warnings import unittest +import importlib +@unittest.skipUnless(importlib.util.find_spec("pyoptsparse") is not None, "pyoptsparse is not installed") class IdleDescentTestCase(unittest.TestCase): def setUp(self): core_subsystems = [prop, geom, aero] From 591a6713da4b4e186cbeb11c4902d12b0d075c7d Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 20 Feb 2024 09:10:37 -0800 Subject: [PATCH 31/34] reverted test values --- .../mission/gasp_based/test/test_idle_descent_estimation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/aviary/mission/gasp_based/test/test_idle_descent_estimation.py b/aviary/mission/gasp_based/test/test_idle_descent_estimation.py index d1363e8c9..c748d6d4b 100644 --- a/aviary/mission/gasp_based/test/test_idle_descent_estimation.py +++ b/aviary/mission/gasp_based/test/test_idle_descent_estimation.py @@ -29,8 +29,8 @@ def test_case1(self): results = descent_range_and_fuel(ode_args=ode_args)['refined_guess'] # Values obtained by running idle_descent_estimation - assert_near_equal(results['distance_flown'], 91.50393932590772, tol) - assert_near_equal(results['fuel_burned'], 235.39517919329228, tol) + assert_near_equal(results['distance_flown'], 91.8911599691433, tol) + assert_near_equal(results['fuel_burned'], 236.73893823639082, tol) if __name__ == "__main__": From 7b09c6ec18dc3aa004bd786eaf6c357dfd96da5e Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 20 Feb 2024 11:51:12 -0800 Subject: [PATCH 32/34] minor changes to tests --- .../test/test_time_integration_phases.py | 26 +++++++++---------- 1 file changed, 12 insertions(+), 14 deletions(-) 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 index 376db30fa..4665096cf 100644 --- a/aviary/mission/flops_based/phases/test/test_time_integration_phases.py +++ b/aviary/mission/flops_based/phases/test/test_time_integration_phases.py @@ -11,7 +11,7 @@ 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 phase_info, aero, prop, geom +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 @@ -23,10 +23,8 @@ @unittest.skipUnless(importlib.util.find_spec("pyoptsparse") is not None, "pyoptsparse is not installed") -class IdleDescentTestCase(unittest.TestCase): +class HE_SGMDescentTestCase(unittest.TestCase): def setUp(self): - core_subsystems = [prop, geom, aero] - 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") @@ -35,9 +33,8 @@ def setUp(self): 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=core_subsystems) - engine = EngineDeck(options=aviary_inputs) - preprocess_propulsion(aviary_inputs, [engine]) + 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'}} @@ -137,25 +134,26 @@ def run_simulation(self, phases, initial_values: dict): # "traj.altitude_initial": {'val': 0, 'units': "ft"}, # "traj.mass_initial": {'val': 171000, 'units': "lbm"}, # "traj.distance_initial": {'val': 0, 'units': "NM"}, - # "traj.velocity": {'val': 0, 'units': "m/s"}, + # "traj.velocity": {'val': .1, 'units': "m/s"}, # } # ode_args = self.ode_args # ode_args['friction_key'] = Mission.Takeoff.ROLLING_FRICTION_COEFFICIENT - # increased_alt_Takeoff = SGMDetailedTakeoff( + # brake_release_to_decision = SGMDetailedTakeoff( # ode_args, - # simupy_args=dict(verbosity=Verbosity.QUIET,) + # simupy_args=dict(verbosity=Verbosity.DEBUG,) # ) - # increased_alt_Takeoff.clear_triggers() - # increased_alt_Takeoff.add_trigger(Dynamic.Mission.ALTITUDE, 500, units='ft') + # brake_release_to_decision.clear_triggers() + # brake_release_to_decision.add_trigger(Dynamic.Mission.VELOCITY, value=167.85, units='kn') # phases = {'HE': { - # 'ode': increased_alt_Takeoff, + # '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['altitude'], 500, self.tol) + # assert_near_equal(final_states['velocity'], 167.85, self.tol) def test_cruise(self): initial_values_cruise = { From a4f04ed607f55b85de15d9858bd432c281093776 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 20 Feb 2024 12:02:10 -0800 Subject: [PATCH 33/34] fixed value when clearing triggers --- aviary/mission/gasp_based/ode/time_integration_base_classes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 47ae1ec5b..4856e19f4 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -395,7 +395,7 @@ def add_trigger(self, state, value, units=None, channel_name=None): def clear_triggers(self): self.triggers = [] self.event_channel_names = [] - self.num_events = [] + self.num_events = 0 def event_equation_function(self, t, x): self.output_equation_function(t, x) From d5a24ebf73cb0c26614b83d35e57c9427a7de4b6 Mon Sep 17 00:00:00 2001 From: Carl Recine Date: Tue, 20 Feb 2024 14:07:04 -0800 Subject: [PATCH 34/34] PR feedback --- aviary/mission/flops_based/ode/mission_ODE.py | 24 --------- .../ode/time_integration_base_classes.py | 51 ++++++++++--------- .../benchmark_tests/test_bench_GwGm.py | 2 +- 3 files changed, 27 insertions(+), 50 deletions(-) diff --git a/aviary/mission/flops_based/ode/mission_ODE.py b/aviary/mission/flops_based/ode/mission_ODE.py index 7b97787ea..87756a921 100644 --- a/aviary/mission/flops_based/ode/mission_ODE.py +++ b/aviary/mission/flops_based/ode/mission_ODE.py @@ -225,30 +225,6 @@ def setup(self): print_level = 0 if analysis_scheme is AnalysisScheme.SHOOTING else 2 - if analysis_scheme is AnalysisScheme.SHOOTING and False: - from aviary.utils.functions import create_printcomp - dummy_comp = create_printcomp( - all_inputs=[ - 't_curr', - Mission.Design.RESERVE_FUEL, - Dynamic.Mission.MASS, - Dynamic.Mission.DISTANCE, - Dynamic.Mission.ALTITUDE, - Dynamic.Mission.FLIGHT_PATH_ANGLE, - ], - input_units={ - 't_curr': 's', - Dynamic.Mission.FLIGHT_PATH_ANGLE: 'deg', - Dynamic.Mission.DISTANCE: 'NM', - }) - self.add_subsystem( - "dummy_comp", - dummy_comp(), - promotes_inputs=["*"],) - self.set_input_defaults( - Dynamic.Mission.DISTANCE, val=0, units='NM') - self.set_input_defaults('t_curr', val=0, units='s') - self.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, atol=1.0e-10, rtol=1.0e-10, 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 4856e19f4..b0e823843 100644 --- a/aviary/mission/gasp_based/ode/time_integration_base_classes.py +++ b/aviary/mission/gasp_based/ode/time_integration_base_classes.py @@ -218,6 +218,7 @@ def __init__( self.t_name = t_name self.states = states + self.state_names = list(states.keys()) self.parameters = parameters self.outputs = outputs @@ -578,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: @@ -641,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: @@ -698,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, - list(next_prob.states.keys()), + next_prob.state_names, return_format='array' ).squeeze() @@ -745,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], - list(prob.states.keys()), + prob.state_names, return_format='array' ) @@ -786,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 @@ -811,7 +812,7 @@ def compute_partials(self, inputs, J): dh_j_dx = prob.compute_totals( [state_name], - list(prob.states.keys()), + prob.state_names, return_format='array').squeeze() dh_dparam[state_idx, :] = prob.compute_totals( @@ -820,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[ @@ -842,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, - list(prob.states.keys()), + prob.state_names, return_format='array').T if param_dict: df_dparam_data[idx, ...] = prob.compute_totals( @@ -957,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...", @@ -966,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], @@ -1063,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/validation_cases/benchmark_tests/test_bench_GwGm.py b/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py index 6ccf14ecf..5e6c5191e 100644 --- a/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py +++ b/aviary/validation_cases/benchmark_tests/test_bench_GwGm.py @@ -61,7 +61,7 @@ def test_bench_GwGm_SNOPT(self): assert_near_equal(prob.get_val("traj.desc2.timeseries.distance")[-1], 3675.0, tolerance=rtol) - @require_pyoptsparse(optimizer="SNOPT") + @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',