Aviary Model Without AviaryProblem#
Because Aviary is, at its core, a very sophisticated wrapper for creating an OpenMDAO problem, an advanced user can replicate Aviary’s capabilities in a custom Python script. At this level, users have full access to Python and OpenMDAO methods that Aviary calls. They can use the complete set of Aviary’s methods, functionalities, and classes to construct and fine-tune their aircraft models.
This enables users to have supreme control over every aspect of the model, including subsystems, connections, and advanced optimization techniques.
We will show this approach using an example that implements all methods that are available in aviary/core/aviary_problem.py but without directly calling AviaryProblem.
Note
For each of these examples we have set max_iter = 0, which means that the optimization will not run. This is done to reduce the computational time for the examples. If you want to run the optimization, you can set max_iter = 100 or some similar value.
# Testing Cell
import os
from aviary.utils.functions import get_path
folder = 'validation_cases/benchmark_tests'
file = 'test_FLOPS_based_sizing_N3CC.py'
get_path(os.path.join(folder, file))
PosixPath('/home/runner/work/Aviary/Aviary/aviary/validation_cases/benchmark_tests/test_FLOPS_based_sizing_N3CC.py')
Advanced Single Aisle Example#
Earlier, we have shown how to follow the standard steps to build (and slightly customize) an Aviary problem, but sometimes you may find that the model you have in mind does not match that predefined structure. If your needs differ from the setup Aviary uses, we can write special methods ourselves. This example will show you how to do that.
In aviary/validation_cases/benchmark_tests folder, there is an aircraft full mission benchmark test test_FLOPS_based_sizing_N3CC.py. Now, we will show how to recreate that example, without using the AviaryProblem. The key is that we follow the same steps:
__init__()load_inputs()check_and_preprocess_inputs()build_model()add_pre_mission_systems()add_phases()add_post_mission_systems()link_phases()add_driver()add_design_variables()(optional)add_objective()setup()run_aviary_problem()
Examining the design case#
"""
This is an example of running a coupled aircraft design-mission optimization in Aviary without using
the lower level APIs.
It runs the same aircraft and mission as the `simple mission` example.
The aircraft is loaded from .csv and the default energy-state phase_info dictionary is imported
from the file.
This script unwraps the subsystem and trajectory builders, exposing how Aviary interacts with
openMDAO and Dymos. It is divided into sections separated by '####' with the level2 api calls
commented out so you can see what happens in each level2 method.
"""
import warnings
import dymos as dm
import openmdao.api as om
import aviary.api as av
from aviary.core.pre_mission_group import PreMissionGroup
from aviary.mission.flight_phase_builder import FlightPhaseOptions
from aviary.mission.energy_state.ode.energy_state_ODE import EnergyStateODE
from aviary.models.missions.energy_state_default import phase_info
from aviary.utils.aviary_values import AviaryValues
from aviary.variable_info.enums import Verbosity
from aviary.variable_info.functions import setup_model_options, setup_trajectory_params
from aviary.variable_info.variable_meta_data import _MetaData as BaseMetaData
from aviary.variable_info.variables import Aircraft, Dynamic, Mission
class L3SubsystemsGroup(om.Group):
"""Group that contains all pre-mission groups of core Aviary subsystems (geometry, mass, propulsion, aerodynamics)."""
def initialize(self):
self.options.declare(
'aviary_options',
types=AviaryValues,
desc='collection of Aircraft/Mission specific options',
)
self.code_origin_overrides = []
# Toggle this boolean option to run with shooting vs collocation transcription:
# shooting = True
shooting = False
prob = av.AviaryProblem()
#####
# prob.load_inputs(csv_path, phase_info)
csv_path = 'models/aircraft/advanced_single_aisle/advanced_single_aisle_FLOPS.csv'
aviary_inputs, _ = av.create_vehicle(csv_path)
engine = av.build_engine_deck(aviary_inputs)
prob.model.mission_info = {}
for phase_name in phase_info:
if phase_name not in ['pre_mission', 'post_mission']:
prob.model.mission_info[phase_name] = phase_info[phase_name]
aviary_inputs.set_val(Mission.Summary.RANGE, 1906.0, units='NM')
prob.require_range_residual = True
prob.target_range = 1906.0
#####
# prob.check_and_preprocess_inputs()
av.preprocess_options(aviary_inputs, engine_models=[engine])
#####
# prob.add_pre_mission_systems()
aerodynamics = av.CoreAerodynamicsBuilder(code_origin=av.LegacyCode('FLOPS'))
geometry = av.CoreGeometryBuilder(code_origin=av.LegacyCode('FLOPS'))
mass = av.CoreMassBuilder(code_origin=av.LegacyCode('FLOPS'))
propulsion = av.CorePropulsionBuilder(engine_models=engine)
prob.model.core_subsystems = {
'propulsion': propulsion,
'geometry': geometry,
'mass': mass,
'aerodynamics': aerodynamics,
}
prob.meta_data = BaseMetaData.copy()
prob.model.add_subsystem(
'pre_mission',
PreMissionGroup(),
promotes_inputs=['aircraft:*', 'mission:*'],
promotes_outputs=['aircraft:*', 'mission:*'],
)
# This is a combination of prob.add_pre_mission_systems and prob.setup()
# In the aviary code add_pre_mission_systems only instantiates the objects and methods, the build method is called in prob.setup()
prob.model.pre_mission.add_subsystem(
'core_propulsion',
propulsion.build_pre_mission(aviary_inputs),
)
# adding another group subsystem to match the L2 example
prob.model.pre_mission.add_subsystem(
'core_subsystems',
L3SubsystemsGroup(aviary_options=aviary_inputs),
promotes_inputs=['*'],
promotes_outputs=['*'],
)
prob.model.pre_mission.core_subsystems.add_subsystem(
'core_geometry',
geometry.build_pre_mission(aviary_inputs),
promotes_inputs=['*'],
promotes_outputs=['*'],
)
prob.model.pre_mission.core_subsystems.add_subsystem(
'core_aerodynamics',
aerodynamics.build_pre_mission(aviary_inputs),
promotes_inputs=['*'],
promotes_outputs=['*'],
)
prob.model.pre_mission.core_subsystems.add_subsystem(
'core_mass',
mass.build_pre_mission(aviary_inputs),
promotes_inputs=['*'],
promotes_outputs=['*'],
)
#####
# prob.add_phases()
phase_list = ['climb', 'cruise', 'descent']
prob.traj = prob.model.add_subsystem('traj', dm.Trajectory())
default_mission_subsystems = [
prob.model.core_subsystems['aerodynamics'],
prob.model.core_subsystems['propulsion'],
]
phases = {}
for phase_idx, phase_name in enumerate(phase_list):
base_phase_options = prob.model.mission_info[phase_name]
phase_options = {}
for key, val in base_phase_options.items():
phase_options[key] = val
phase_options['user_options'] = {}
for key, val in base_phase_options['user_options'].items():
phase_options['user_options'][key] = val
# Create the dymos phase object and add it to the trajectory
# phase_builder = EnergyPhase
# Unpack the phase info without using the phase builder
# phase_object = phase_builder.from_phase_info(
# phase_name, phase_options, default_mission_subsystems, meta_data=prob.meta_data
# )
# user_options dict entries need reformatting
# if the value is not a tuple, wrap it in a tuple with the second entry of 'unitless'
for key, value in phase_options['user_options'].items():
if not isinstance(value, tuple):
phase_options['user_options'][key] = (value, 'unitless')
subsystem_options = phase_options.get('subsystem_options', {})
user_options = phase_options.get('user_options', ())
initial_guesses = AviaryValues(
phase_options.get('initial_guesses', ())
) # Not used in this example
# instantiate the PhaseBuilderBaseClass:
# phase_builder = cls(
# name,
# subsystem_options=subsystem_options,
# user_options=user_options,
# initial_guesses=initial_guesses,
# meta_data=meta_data,
# core_subsystems=core_subsystems,
# external_subsystems=external_subsystems,
# transcription=transcription,
# )
# This basically just adds these objects to the class - we have them all available in this L3 script so have no need for the extra class!
transcription = None
ode_class = None
is_analytic_phase = False
num_nodes = 5 # this is only used for analytic phases
# Now build the phase using the instantiated phase object:
# phase = phase_object.build_phase(aviary_options=aviary_inputs)
# phase: dm.Phase = super().build_phase(aviary_options)
# For L3 we need to do this without the builder
ode_class = EnergyStateODE
user_options = FlightPhaseOptions(user_options)
num_segments = user_options['num_segments']
order = user_options['order']
if shooting:
nodes_per_seg = order * num_segments
transcription = dm.PicardShooting(
num_segments=1, nodes_per_seg=nodes_per_seg, solve_segments='forward'
)
else:
seg_ends, _ = dm.utils.lgl.lgl(num_segments + 1)
transcription = dm.Radau(
num_segments=num_segments, order=order, compressed=True, segment_ends=seg_ends
)
kwargs = {
'meta_data': prob.meta_data,
'subsystem_options': subsystem_options,
'throttle_enforcement': user_options['throttle_enforcement'],
'throttle_allocation': user_options['throttle_allocation'],
'subsystems': default_mission_subsystems,
}
kwargs = {'aviary_options': aviary_inputs, **kwargs}
phase = dm.Phase(ode_class=ode_class, transcription=transcription, ode_init_kwargs=kwargs)
tx_mission_bus = dm.GaussLobatto(
num_segments=transcription.options['num_segments'], order=3, compressed=True
) # equally spaced points in this output - this is going to be better than using picard or radau.
phase.add_timeseries(name='mission_bus_variables', transcription=tx_mission_bus, subset='all')
num_engine_type = len(
aviary_inputs.get_val(Aircraft.Engine.NUM_ENGINES)
) # not used in this example
throttle_enforcement = user_options['throttle_enforcement']
no_descent = user_options['no_descent'] # not used in this example
no_climb = user_options['no_climb'] # not used in this example
constraints = user_options['constraints'] # not used in this example
ground_roll = user_options['ground_roll'] # not used in this example
phase.add_state(
name=Dynamic.Vehicle.MASS,
fix_initial=False,
input_initial=False,
lower=user_options['mass_bounds'][0][0],
upper=user_options['mass_bounds'][0][1],
units=user_options['mass_bounds'][1],
rate_source=Dynamic.Vehicle.Propulsion.FUEL_FLOW_RATE_NEGATIVE_TOTAL,
ref=user_options['mass_ref'][0],
ref0=user_options['mass_ref0'][0],
defect_ref=user_options['mass_defect_ref'][0],
solve_segments='forward' if user_options['mass_solve_segments'] else None,
)
phase.add_state(
name=Dynamic.Mission.DISTANCE,
fix_initial=False,
input_initial=False,
lower=user_options['distance_bounds'][0][0],
upper=user_options['distance_bounds'][0][1],
units=user_options['distance_bounds'][1],
rate_source=Dynamic.Mission.DISTANCE_RATE,
ref=user_options['distance_ref'][0],
ref0=user_options['distance_ref0'][0],
defect_ref=user_options['distance_defect_ref'][0],
solve_segments='forward' if user_options['distance_solve_segments'] else None,
)
phase.add_control(
Dynamic.Atmosphere.MACH,
targets=Dynamic.Atmosphere.MACH,
opt=user_options[f'mach_optimize'],
lower=user_options[f'mach_bounds'][0][0],
upper=user_options[f'mach_bounds'][0][1],
ref=user_options['mach_ref'][0],
ref0=user_options['mach_ref0'][0],
control_type='polynomial',
order=user_options[f'mach_polynomial_order'],
rate_targets=[Dynamic.Atmosphere.MACH_RATE],
)
phase.add_timeseries_output(Dynamic.Atmosphere.MACH)
phase.add_control(
Dynamic.Mission.ALTITUDE,
targets=Dynamic.Mission.ALTITUDE,
opt=user_options[f'altitude_optimize'],
lower=user_options[f'altitude_bounds'][0][0],
upper=user_options[f'altitude_bounds'][0][1],
ref=user_options['altitude_ref'][0],
ref0=user_options['altitude_ref0'][0],
control_type='polynomial',
order=user_options[f'altitude_polynomial_order'],
units=user_options[f'altitude_bounds'][1],
rate_targets=[Dynamic.Mission.ALTITUDE_RATE],
)
phase.add_timeseries_output(Dynamic.Mission.ALTITUDE)
# if throttle_enforcement == 'control':
# add throttle as control
# not used in this example
phase.add_timeseries_output(
Dynamic.Vehicle.Propulsion.THRUST_TOTAL,
output_name=Dynamic.Vehicle.Propulsion.THRUST_TOTAL,
units='lbf',
)
phase.add_timeseries_output(Dynamic.Vehicle.DRAG, output_name=Dynamic.Vehicle.DRAG, units='lbf')
phase.add_timeseries_output(
Dynamic.Mission.SPECIFIC_ENERGY_RATE_EXCESS,
output_name=Dynamic.Mission.SPECIFIC_ENERGY_RATE_EXCESS,
units='m/s',
)
phase.add_timeseries_output(
Dynamic.Vehicle.Propulsion.FUEL_FLOW_RATE_NEGATIVE_TOTAL,
output_name=Dynamic.Vehicle.Propulsion.FUEL_FLOW_RATE_NEGATIVE_TOTAL,
units='lbm/h',
)
phase.add_timeseries_output(
Dynamic.Vehicle.Propulsion.ELECTRIC_POWER_IN_TOTAL,
output_name=Dynamic.Vehicle.Propulsion.ELECTRIC_POWER_IN_TOTAL,
units='kW',
)
phase.add_timeseries_output(
Dynamic.Mission.ALTITUDE_RATE,
output_name=Dynamic.Mission.ALTITUDE_RATE,
units='ft/s',
)
if throttle_enforcement != 'control':
phase.add_timeseries_output(
Dynamic.Vehicle.Propulsion.THROTTLE,
output_name=Dynamic.Vehicle.Propulsion.THROTTLE,
units='unitless',
)
phase.add_timeseries_output(
Dynamic.Mission.VELOCITY,
output_name=Dynamic.Mission.VELOCITY,
units='m/s',
)
phase.add_path_constraint(
Dynamic.Vehicle.Propulsion.THROTTLE,
lower=0.0,
upper=1.0,
units='unitless',
)
prob.traj.add_phase(phase_name, phase)
phases[phase_name] = phase
climb = phases['climb']
cruise = phases['cruise']
descent = phases['descent']
externs = {'climb': {}, 'cruise': {}, 'descent': {}}
for default_subsys in default_mission_subsystems:
params = default_subsys.get_parameters(aviary_inputs=aviary_inputs, phase_info={})
for key, val in params.items():
for phname in externs:
externs[phname][key] = val
prob.traj = setup_trajectory_params(
prob.model, prob.traj, aviary_inputs, external_parameters=externs
)
# Need aviary inputs assigned to the problem object for other functions below
prob.aviary_inputs = aviary_inputs
#####
# prob.add_post_mission_systems()
prob.model.add_subsystem(
'post_mission',
om.Group(),
promotes_inputs=['*'],
promotes_outputs=['*'],
)
climb.set_state_options(Dynamic.Vehicle.MASS, fix_initial=False, input_initial=False)
climb.set_state_options(Dynamic.Mission.DISTANCE, fix_initial=True, input_initial=False)
climb.set_time_options(
fix_initial=False,
initial_bounds=(0, 0),
initial_ref=600,
duration_bounds=(3840, 11520),
duration_ref=7680.0,
)
cruise.set_time_options(
duration_bounds=(3390, 18000),
duration_ref=10695.0,
)
descent.set_time_options(
duration_bounds=(1740, 5220),
duration_ref=3480.0,
)
eq = prob.model.add_subsystem(
f'link_climb_mass',
om.EQConstraintComp(),
promotes_inputs=[('rhs:mass', Mission.Summary.GROSS_MASS)],
)
eq.add_eq_output('mass', eq_units='lbm', normalize=False, ref=100000.0, add_constraint=True)
prob.model.connect(
f'traj.climb.states:mass',
f'link_climb_mass.lhs:mass',
src_indices=[0],
flat_src_indices=True,
)
prob.model.add_subsystem(
'range_constraint',
om.ExecComp(
'range_resid = target_range - actual_range',
target_range={'val': prob.target_range, 'units': 'NM'},
actual_range={'val': prob.target_range, 'units': 'NM'},
range_resid={'val': 30, 'units': 'NM'},
),
promotes_inputs=[
('actual_range', Mission.Summary.RANGE),
'target_range',
],
promotes_outputs=[('range_resid', Mission.Constraints.RANGE_RESIDUAL)],
)
prob.model.add_constraint(Mission.Constraints.MASS_RESIDUAL, equals=0.0, ref=1.0e5)
# for reference this is the end of builder.add_post_mission_systems()
ecomp = om.ExecComp(
'fuel_burned = initial_mass - mass_final',
initial_mass={'units': 'lbm'},
mass_final={'units': 'lbm'},
fuel_burned={'units': 'lbm'},
)
prob.model.post_mission.add_subsystem(
'fuel_burned',
ecomp,
promotes=[('fuel_burned', Mission.Summary.FUEL_BURNED)],
)
prob.model.connect(
f'traj.climb.timeseries.mass',
'fuel_burned.initial_mass',
src_indices=[0],
)
prob.model.connect(
f'traj.descent.timeseries.mass',
'fuel_burned.mass_final',
src_indices=[-1],
)
reserve_fuel_additional = prob.aviary_inputs.get_val(
Aircraft.Design.RESERVE_FUEL_ADDITIONAL, units='lbm'
)
reserve_fuel = om.ExecComp(
'reserve_fuel = reserve_fuel_frac_mass + reserve_fuel_additional + reserve_fuel_burned',
reserve_fuel={'units': 'lbm', 'shape': 1},
reserve_fuel_frac_mass={'units': 'lbm', 'val': 0},
reserve_fuel_additional={'units': 'lbm', 'val': reserve_fuel_additional},
reserve_fuel_burned={'units': 'lbm', 'val': 0},
)
prob.model.post_mission.add_subsystem(
'reserve_fuel',
reserve_fuel,
promotes_inputs=[
'reserve_fuel_frac_mass',
('reserve_fuel_additional', Aircraft.Design.RESERVE_FUEL_ADDITIONAL),
('reserve_fuel_burned', Mission.Summary.RESERVE_FUEL_BURNED),
],
promotes_outputs=[('reserve_fuel', Mission.Design.RESERVE_FUEL)],
)
ecomp = om.ExecComp(
'overall_fuel = fuel_burned + reserve_fuel',
overall_fuel={'units': 'lbm', 'shape': 1},
fuel_burned={'units': 'lbm'}, # from regular_phases only
reserve_fuel={'units': 'lbm', 'shape': 1},
)
prob.model.post_mission.add_subsystem(
'fuel_calc',
ecomp,
promotes_inputs=[
('fuel_burned', Mission.Summary.FUEL_BURNED),
('reserve_fuel', Mission.Design.RESERVE_FUEL),
],
promotes_outputs=[('overall_fuel', Mission.Summary.TOTAL_FUEL_MASS)],
)
# If target distances have been set per phase then there is a block of code to add here.
# In this case individual phases don't have target distances.
ecomp = om.ExecComp(
'mass_resid = operating_empty_mass + overall_fuel + payload_mass - initial_mass',
operating_empty_mass={'units': 'lbm'},
overall_fuel={'units': 'lbm'},
payload_mass={'units': 'lbm'},
initial_mass={'units': 'lbm'},
mass_resid={'units': 'lbm'},
)
prob.model.post_mission.add_subsystem(
'mass_constraint',
ecomp,
promotes_inputs=[
('operating_empty_mass', Mission.Summary.OPERATING_MASS),
('overall_fuel', Mission.Summary.TOTAL_FUEL_MASS),
('payload_mass', Aircraft.CrewPayload.TOTAL_PAYLOAD_MASS),
('initial_mass', Mission.Summary.GROSS_MASS),
],
promotes_outputs=[('mass_resid', Mission.Constraints.MASS_RESIDUAL)],
)
ecomp = om.ExecComp(
'excess_fuel_capacity = total_fuel_capacity - unusable_fuel - overall_fuel',
total_fuel_capacity={'units': 'lbm'},
unusable_fuel={'units': 'lbm'},
overall_fuel={'units': 'lbm'},
excess_fuel_capacity={'units': 'lbm'},
)
prob.model.post_mission.add_subsystem(
'excess_fuel_constraint',
ecomp,
promotes_inputs=[
('total_fuel_capacity', Aircraft.Fuel.TOTAL_CAPACITY),
('unusable_fuel', Aircraft.Fuel.UNUSABLE_FUEL_MASS),
('overall_fuel', Mission.Summary.TOTAL_FUEL_MASS),
],
promotes_outputs=[('excess_fuel_capacity', Mission.Constraints.EXCESS_FUEL_CAPACITY)],
)
prob.model.add_constraint(Mission.Constraints.EXCESS_FUEL_CAPACITY, lower=0, units='lbm')
#####
# prob.link_phases()
all_subsystems = []
all_subsystems.append(prob.model.core_subsystems['propulsion'])
prob.traj.link_phases(phase_list, ['time'], ref=None, connected=True)
prob.traj.link_phases(phase_list, [Dynamic.Vehicle.MASS], ref=None, connected=True)
prob.traj.link_phases(phase_list, [Dynamic.Mission.DISTANCE], ref=None, connected=True)
prob.model.connect(
f'traj.descent.timeseries.distance',
Mission.Summary.RANGE,
src_indices=[-1],
flat_src_indices=True,
)
#####
# Each driver requires different settings - uncomment the one to be used:
# prob.add_driver('SLSQP', max_iter=50)
# SLSQP Optimizer Settings
# prob.driver = om.ScipyOptimizeDriver()
# prob.driver.options['optimizer'] = 'SLSQP'
# prob.driver.declare_coloring(show_summary=False)
# prob.driver.options['disp'] = True
# prob.driver.options['tol'] = 1e-9
# prob.driver.options['maxiter'] = 50
# prob.add_driver('IPOPT', max_iter=50)
# IPOPT Optimizer Settings
prob.driver = om.pyOptSparseDriver()
prob.driver.options['optimizer'] = 'IPOPT'
prob.driver.declare_coloring(show_summary=False)
prob.driver.opt_settings['print_user_options'] = 'no'
prob.driver.opt_settings['print_frequency_iter'] = 10
prob.driver.opt_settings['print_level'] = 3
prob.driver.opt_settings['tol'] = 1.0e-6
prob.driver.opt_settings['mu_init'] = 1e-5
prob.driver.opt_settings['max_iter'] = 300
prob.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
prob.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
prob.driver.opt_settings['mu_strategy'] = 'monotone'
prob.driver.options['print_results'] = 'minimal'
# prob.add_driver('SNOPT', max_iter=50)
# SNOPT Optimizer Settings #
# prob.driver = om.pyOptSparseDriver()
# prob.driver.options['optimizer'] = 'SNOPT'
# prob.driver.declare_coloring(show_summary=False)
# prob.driver.opt_settings['iSumm'] = 6
# prob.driver.opt_settings['iPrint'] = 0
# prob.driver.opt_settings['Major iterations limit'] = 50
# prob.driver.opt_settings['Major optimality tolerance'] = 1e-4
# prob.driver.opt_settings['Major feasibility tolerance'] = 1e-7
# prob.driver.options['print_results'] = 'minimal'
#####
# prob.add_design_variables()
prob.model.add_design_var(
Mission.Design.GROSS_MASS,
lower=100000.0,
upper=None,
units='lbm',
ref=175e3,
)
prob.model.add_design_var(
Mission.Summary.GROSS_MASS,
lower=100000.0,
upper=None,
units='lbm',
ref=175e3,
)
prob.model.add_subsystem(
'gtow_constraint',
om.EQConstraintComp(
'GTOW',
eq_units='lbm',
normalize=True,
add_constraint=True,
),
promotes_inputs=[
('lhs:GTOW', Mission.Design.GROSS_MASS),
('rhs:GTOW', Mission.Summary.GROSS_MASS),
],
)
prob.model.add_constraint(Mission.Constraints.RANGE_RESIDUAL, equals=0, ref=10)
#####
# prob.add_objective()
prob.model.add_subsystem(
'fuel_obj',
om.ExecComp(
'reg_objective = overall_fuel/10000 + ascent_duration/30.',
reg_objective={'val': 0.0, 'units': 'unitless'},
ascent_duration={'units': 's', 'shape': 1},
overall_fuel={'units': 'lbm'},
),
promotes_inputs=[
('ascent_duration', Mission.Takeoff.ASCENT_DURATION),
('overall_fuel', Mission.Summary.TOTAL_FUEL_MASS),
],
promotes_outputs=[('reg_objective', Mission.Objectives.FUEL)],
)
prob.model.add_objective(Mission.Objectives.FUEL, ref=1)
prob.model.add_subsystem(
'range_obj',
om.ExecComp(
'reg_objective = -actual_range/1000 + ascent_duration/30.',
reg_objective={'val': 0.0, 'units': 'unitless'},
ascent_duration={'units': 's', 'shape': 1},
actual_range={'val': prob.target_range, 'units': 'NM'},
),
promotes_inputs=[
('actual_range', Mission.Summary.RANGE),
('ascent_duration', Mission.Takeoff.ASCENT_DURATION),
],
promotes_outputs=[('reg_objective', Mission.Objectives.RANGE)],
)
#####
# prob.setup()
setup_model_options(prob, prob.aviary_inputs, prob.meta_data)
with warnings.catch_warnings():
prob.model.aviary_inputs = prob.aviary_inputs
prob.model.meta_data = prob.meta_data
with warnings.catch_warnings():
warnings.simplefilter('ignore', om.OpenMDAOWarning)
warnings.simplefilter('ignore', om.PromotionWarning)
om.Problem.setup(prob, check=False)
# set initial guesses manually
control_keys = ['mach', 'altitude']
state_keys = ['mass', Dynamic.Mission.DISTANCE]
guesses = {}
guesses['mach_climb'] = ([0.2, 0.72], 'unitless')
guesses['altitude_climb'] = ([0, 32000.0], 'ft')
guesses['time_climb'] = ([0, 3840.0], 's')
guesses['mach_cruise'] = ([0.72, 0.72], 'unitless')
guesses['altitude_cruise'] = ([32000.0, 34000.0], 'ft')
guesses['time_cruise'] = ([3840.0, 3390.0], 's')
guesses['mach_descent'] = ([0.72, 0.36], 'unitless')
guesses['altitude_descent'] = ([34000.0, 500.0], 'ft')
guesses['time_descent'] = ([7230.0, 1740.0], 's')
climb.set_time_val(
initial=guesses['time_climb'][0][0], duration=guesses['time_climb'][0][1], units='s'
)
climb.set_control_val('mach', vals=guesses['mach_climb'][0], time_vals=[-1, 1], units='unitless')
climb.set_control_val('altitude', vals=guesses['altitude_climb'][0], time_vals=[-1, 1], units='ft')
climb.set_state_val('mass', 125000, units='lbm')
cruise.set_time_val(
initial=guesses['time_cruise'][0][0], duration=guesses['time_cruise'][0][1], units='s'
)
cruise.set_control_val('mach', vals=guesses['mach_cruise'][0], time_vals=[-1, 1], units='unitless')
cruise.set_control_val(
'altitude', vals=guesses['altitude_cruise'][0], time_vals=[-1, 1], units='ft'
)
cruise.set_state_val('mass', 125000, units='lbm')
descent.set_time_val(
initial=guesses['time_descent'][0][0], duration=guesses['time_descent'][0][1], units='s'
)
descent.set_control_val(
'mach', vals=guesses['mach_descent'][0], time_vals=[-1, 1], units='unitless'
)
descent.set_control_val(
'altitude', vals=guesses['altitude_descent'][0], time_vals=[-1, 1], units='ft'
)
descent.set_state_val('mass', 125000, units='lbm')
prob.set_val(Mission.Design.GROSS_MASS, 175400, units='lbm')
prob.set_val(Mission.Summary.GROSS_MASS, 175400, units='lbm')
prob.verbosity = Verbosity.BRIEF
prob.run_aviary_problem()
# Uncomment these lines to get printouts of every variable in the openmdao model
# prob.model.list_vars(units=True, print_arrays=True)
# prob.list_driver_vars(print_arrays=True)
Warning: Invalid options for non-optimal control 'mach' in phase 'climb': lower, upper, ref
Warning: Invalid options for non-optimal control 'altitude' in phase 'climb': lower, upper, ref
Warning: Invalid options for non-optimal control 'mach' in phase 'cruise': lower, upper, ref
Warning: Invalid options for non-optimal control 'altitude' in phase 'cruise': lower, upper, ref
Warning: Invalid options for non-optimal control 'mach' in phase 'descent': lower, upper, ref
Warning: Invalid options for non-optimal control 'altitude' in phase 'descent': lower, upper, ref
/usr/share/miniconda/envs/test/lib/python3.13/site-packages/openmdao/utils/relevance.py:1296: OpenMDAOWarning:The following groups have a nonlinear solver that computes gradients and will be treated as atomic for the purposes of determining which systems are included in the optimization iteration:
traj.phases.climb.rhs_all.solver_sub
traj.phases.cruise.rhs_all.solver_sub
traj.phases.descent.rhs_all.solver_sub
/usr/share/miniconda/envs/test/lib/python3.13/site-packages/openmdao/solvers/linear/linear_rhs_checker.py:175: SolverWarning:DirectSolver in 'traj.phases.cruise.indep_states' <class StateIndependentsComp>: 'rhs_checking' is active but no redundant adjoint dependencies were found, so caching has been disabled.
/usr/share/miniconda/envs/test/lib/python3.13/site-packages/openmdao/solvers/linear/linear_rhs_checker.py:175: SolverWarning:DirectSolver in 'traj.phases.descent.indep_states' <class StateIndependentsComp>: 'rhs_checking' is active but no redundant adjoint dependencies were found, so caching has been disabled.
/usr/share/miniconda/envs/test/lib/python3.13/site-packages/openmdao/core/total_jac.py:1696: DerivativesWarning:The following design variables have no impact on the constraints or objective at the current design point:
traj.climb.t_initial, inds=[0]
Total number of variables............................: 96
variables with only lower bounds: 93
variables with lower and upper bounds: 3
variables with only upper bounds: 0
Total number of equality constraints.................: 94
Total number of inequality constraints...............: 61
inequality constraints with only lower bounds: 1
inequality constraints with lower and upper bounds: 60
inequality constraints with only upper bounds: 0
Optimization Problem -- Optimization using pyOpt_sparse
================================================================================
Objective Function: _objfunc
Objectives
Index Name Value
0 mission:objectives:fuel 1.323443E+00
Variables (c - continuous, i - integer, d - discrete)
Index Name Type Lower Bound Value Upper Bound Status
2 traj.climb.t_initial_0 c 0.000000E+00 0.000000E+00 0.000000E+00 lu
3 traj.climb.t_duration_0 c 5.000000E-01 5.000006E-01 1.500000E+00 l
Constraints (i - inequality, e - equality)
Index Name Type Lower Value Upper Status Lagrange Multiplier (N/A)
154 traj.descent.throttle[path] i 0.000000E+00 5.221229E-07 1.000000E+00 l 9.00000E+100
Number of Iterations....: 15
(scaled) (unscaled)
Objective...............: 1.3234433662315992e+00 1.3234433662315992e+00
Dual infeasibility......: 3.6142166968691972e-08 3.6142166968691972e-08
Constraint violation....: 2.2737367544323207e-14 2.2737367544323207e-14
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0924048285872868e-08 9.0924048285872868e-08
Overall NLP error.......: 9.0924048285872868e-08 9.0924048285872868e-08
Number of objective function evaluations = 16
Number of objective gradient evaluations = 16
Number of equality constraint evaluations = 16
Number of inequality constraint evaluations = 16
Number of equality constraint Jacobian evaluations = 16
Number of inequality constraint Jacobian evaluations = 16
Number of Lagrangian Hessian evaluations = 0
Total seconds in IPOPT = 5.935
EXIT: Optimal Solution Found.
This model demonstrates the flexibility of a custom Python script. For example, we do not load the aircraft model from a .csv file but from a Python file using the get_flops_inputs() method.
This function not only reads Aviary and mission variables but also builds the engine. More information can be found in models/aircraft/advanced_single_aisle/advanced_single_aisle_data.py.
Note that we can read large single aisle aircraft inputs this way as well:
aviary_inputs = get_flops_inputs('LargeSingleAisle1FLOPS')
aviary_inputs = get_flops_inputs('LargeSingleAisle2FLOPS')
aviary_inputs = get_flops_inputs('LargeSingleAisle2FLOPSdw')
aviary_inputs = get_flops_inputs('LargeSingleAisle2FLOPSalt')
The data files are at:
aviary/models/aircraft/large_single_aisle_1/large_single_aisle_1_FLOPS_data.py
aviary/models/aircraft/large_single_aisle_2/large_single_aisle_2_FLOPS_data.py
aviary/models/aircraft/large_single_aisle_2/large_single_aisle_2_detailwing_FLOPS_data.py
aviary/models/aircraft/large_single_aisle_2/large_single_aisle_2_altwt_FLOPS_data.py
respectively.
Discussing this example in more detail#
We move all the code blocks on taxi to add_pre_mission_systems() function because it is how it is done in aviary/core/aviary_problem.py. Similarly, all the code blocks on landing are moved to add_post_mission_systems() function. Be careful! Generally speaking, not all components can be moved around due to the expected order of execution.
In aviary/validation_cases/benchmark_tests folder, there is another N3CC model test_FLOPS_based_sizing_N3CC.py. If we had started from that model, you would need to have an add_design_variables() function.