Using the Integrator
In this tutorial, we will build off the previous minimal example tutorial by adding a numerical integrator to compute fuel burn throughout the mission. This will require the following additions to the previous aircraft model:
A component to compute fuel flow rate using thrust and thrust-specific fuel consumption (TSFC)
An integrator to integrate the fuel flow rate w.r.t. time to compute the weight of fuel burned at each numerical integration point
A component to compute the weight at each integration point by subtracting the fuel burned from the takeoff weight
Other than these changes to the aircraft model the code will look very similar to the minimal example, so we will gloss over some details that are described more in the minimal example. If you have not done so already, it is recommended to go through that tutorial first.
Note
The script described in this tutorial is called minimal_integrator.py and can be found in the examples folder.
Imports
import openmdao.api as om
from openconcept.examples.minimal import Aircraft, setup_problem # build off this aircraft model
from openconcept.mission import BasicMission
from openconcept.utilities import Integrator
The notable addition to the imports is OpenConcept’s Integrator class. We also import the Aircraft class and setup_problem function from the previous tutorial.
Aircraft model
The aircraft model is no longer a set of explicit equations; it now requires a combination of OpenMDAO components to compute the weight.
For this reason, the aircraft model is now an OpenMDAO Group instead of an ExplicitComponent.
Let’s first take a look at the code for the entire aircraft model and then we will break down what is happening in each section.
class AircraftWithFuelBurn(om.Group):
"""
This model takes the simplified aircraft model from the minimal example, but adds
a fuel flow computation using TSFC and integrates it to compute fuel burn.
"""
# rst Options
def initialize(self):
self.options.declare("num_nodes", default=1, desc="Number of analysis points per phase")
self.options.declare("flight_phase", default=None) # required by OpenConcept but unused in this example
# rst Setup
def setup(self):
nn = self.options["num_nodes"]
# Add the aircraft model from the minimal example to build off of.
# Don't promote the weight from this model because we're going to compute a new
# one using the fuel burn.
# rst Simple aircraft (beg)
self.add_subsystem(
"simple_aircraft",
Aircraft(num_nodes=nn),
promotes_inputs=[
"fltcond|CL",
"throttle",
"fltcond|q",
"ac|geom|wing|S_ref",
"ac|weights|TOW",
"ac|propulsion|max_thrust",
"ac|aero|L_over_D",
],
promotes_outputs=["thrust", "drag"],
)
# rst Simple aircraft (end)
# Use an OpenMDAO ExecComp to compute the fuel flow rate using the thrust and TSFC
# rst Fuel flow (beg)
self.add_subsystem(
"fuel_flow_calc",
om.ExecComp(
"fuel_flow = TSFC * thrust",
fuel_flow={"units": "kg/s", "shape": nn},
TSFC={"units": "kg/N/s", "shape": 1},
thrust={"units": "N", "shape": nn},
has_diag_partials=True,
),
promotes_inputs=[("TSFC", "ac|propulsion|TSFC"), "thrust"],
)
# rst Fuel flow (end)
# Integrate the fuel flow rate to compute fuel burn
# rst Integrator (beg)
integ = self.add_subsystem(
"fuel_integrator", Integrator(num_nodes=nn, diff_units="s", time_setup="duration", method="simpson")
)
integ.add_integrand("fuel_burned", rate_name="fuel_flow", units="kg")
self.connect("fuel_flow_calc.fuel_flow", "fuel_integrator.fuel_flow")
# rst Integrator (end)
# Compute the current weight by subtracting the fuel burned from the takeoff weight
# rst Weight (beg)
self.add_subsystem(
"weight_calc",
om.ExecComp(
"weight = TOW - fuel_burned",
units="kg",
weight={"shape": nn},
TOW={"shape": 1},
fuel_burned={"shape": nn},
has_diag_partials=True,
),
promotes_inputs=[("TOW", "ac|weights|TOW")],
promotes_outputs=["weight"],
)
self.connect("fuel_integrator.fuel_burned", "weight_calc.fuel_burned")
# rst Weight (end)
Options
The options are the same as the minimal example tutorial.
def initialize(self):
self.options.declare("num_nodes", default=1, desc="Number of analysis points per phase")
self.options.declare("flight_phase", default=None) # required by OpenConcept but unused in this example
Setup
Note
The order you add components to OpenMDAO groups (using add_subsystem) matters!
Generally, it is best to try to add components in the order that achieves as much feed-forward variable passing as possible.
For example, we have a component that computes thrust and another that takes thrust as an input.
To make this feed-forward, we add the component that takes thrust as an input after the component that computes it.
Thrust and drag from minimal aircraft
The first step in setting up the new aircraft model is to add the simplified aircraft to the group. We still use this model to compute thrust and drag, but the weight calculation will be modified. For this reason, we promote only the thrust and drag outputs to the new aircraft model group level. All the inputs are still required, so we promote them all to the group level (this way OpenConcept will automatically connect them as we discussed last time). If you are confused about the promotion, check out the OpenMDAO documentation.
self.add_subsystem(
"simple_aircraft",
Aircraft(num_nodes=nn),
promotes_inputs=[
"fltcond|CL",
"throttle",
"fltcond|q",
"ac|geom|wing|S_ref",
"ac|weights|TOW",
"ac|propulsion|max_thrust",
"ac|aero|L_over_D",
],
promotes_outputs=["thrust", "drag"],
)
Fuel flow rate
Next, we need to compute the fuel flow rate to pass to the integrator.
Since this is a simple function of thrust and TSFC, we use an OpenMDAO ExecComp (the OpenMDAO docs are very thorough if you are confused about the syntax).
We give it the fuel flow equation we want it to evaluate and define the units and shape of each parameter.
Notice that fuel flow and thrust are both vectors because they are evaluated at each numerical integration point and will change throughout each flight phase.
The TSFC is a scalar because it is a single constant parameter defined for the aircraft.
Finally, we promote the inputs.
Thrust is automatically connected to the thrust output from the minimal aircraft model.
TSFC is promoted to a name beginning with "ac|" so that the mission analysis promotes the variable to the top level so we can set it the same way as the other aircraft parameters.
self.add_subsystem(
"fuel_flow_calc",
om.ExecComp(
"fuel_flow = TSFC * thrust",
fuel_flow={"units": "kg/s", "shape": nn},
TSFC={"units": "kg/N/s", "shape": 1},
thrust={"units": "N", "shape": nn},
has_diag_partials=True,
),
promotes_inputs=[("TSFC", "ac|propulsion|TSFC"), "thrust"],
)
Integrator
Now we are ready to add the integration.
This is done by adding an OpenConcept Integrator component to the model.
After adding the integrator, we add an integrated variable and associated variable to integrate using the integrator’s add_integrand method.
Let’s step through all the details of these calls—there’s a lot to unpack.
integ = self.add_subsystem(
"fuel_integrator", Integrator(num_nodes=nn, diff_units="s", time_setup="duration", method="simpson")
)
integ.add_integrand("fuel_burned", rate_name="fuel_flow", units="kg")
self.connect("fuel_flow_calc.fuel_flow", "fuel_integrator.fuel_flow")
When Integrator is initialized, there are a few important options that must be set.
As we’ve seen before, we set num_nodes to tell it how many integration points to use.
diff_units are the units of the differential.
For example, in our equation we are computing
The differential is \(dt\) and has units of time (we’ll use seconds here).
The time_setup option sets what information the integrator uses to figure out the time at each integration point.
The options are "dt", "duration", or "bounds".
"dt"creates an input called"dt"that specifies the time spacing between each numerical integration point"duration"creates an input called"duration"that specifies the total time of the phase. The time between each integration point is computed by dividing the duration by the number of time steps (number of nodes minus one). This is the most common choice for the time setup and has the advantage that OpenConcept automatically connects the"duration"input to the mission-level duration, so there is no manual time connection needed."bounds"creates inputs called"t_initial"and"t_final"that specify the initial and final time of the phase. This internally computes duration and then time is computed the same was as for the duration approach.
The final option is the integration scheme.
The two options are "bdf3" and "simpson".
"bdf3" uses the third-order-accurate BDF3 integration scheme.
"simpson" uses Simpson’s rule.
Simpson’s rule is the most common choice for use in OpenConcept.
In the next line we add information about the quantity we want to integrate.
We first define the name of the integrated quantity: "fuel_burned".
This will become a vector output of the integrator (accessed in this case as "fuel_integrator.fuel_burned").
We then define the rate we want integrated: "fuel_flow".
This will create a vector input to the integrator called "fuel_flow".
It also automatically adds an input called "fuel_flow_initial" and an output called "fuel_flow_final".
Instead of appending "_initial" or "_final", these names can be set manually using the start_name and end_name optional arguments.
The final value variable of each phase is automatically linked to the initial value variable of the following one.
The initial value in the first mission phase is zero by default but can be changed either using the start_val optional argument or by setting the variable in a usual OpenMDAO way with an IndepVarComp or set_input_defaults.
We set the units of the fuel burn (the integrated quantity) to kilograms.
Other available options can be found in the Integrator source docs.
The final step is to connect the fuel flow output from the fuel flow computation component to the integrator’s fuel flow input.
Mission
The rest of the code will look very similar to the minimal example.
class MissionAnalysisWithFuelBurn(om.Group):
"""
OpenMDAO group for basic three-phase climb, cruise, descent mission.
The only top-level aircraft variable that the aircraft model uses is
the wing area, so that must be defined.
"""
def initialize(self):
self.options.declare("num_nodes", default=1, desc="Number of analysis points per phase")
def setup(self):
iv = self.add_subsystem("ac_vars", om.IndepVarComp(), promotes_outputs=["*"])
iv.add_output("ac|geom|wing|S_ref", val=25.0, units="m**2")
iv.add_output("ac|weights|TOW", val=5e3, units="kg")
iv.add_output("ac|propulsion|max_thrust", val=1e4, units="N")
iv.add_output("ac|propulsion|TSFC", val=20.0, units="g/kN/s")
iv.add_output("ac|aero|L_over_D", val=10.0)
# Define the mission
self.add_subsystem(
"mission",
BasicMission(aircraft_model=AircraftWithFuelBurn, num_nodes=self.options["num_nodes"]),
promotes_inputs=["ac|*"],
)
The mission is identical except for two changes.
Firstly, we set the TSFC variable called "ac|propulsion|TSFC".
Secondly, the aircraft_model passed to the mission analysis component is now AircraftWithFuelBurn.
Run script
We reuse the setup_problem function from the minimal example.
The remaining code is the same, except for adding a couple more variables of interest to the output plot.
if __name__ == "__main__":
# Process command line argument to optionally not show figures and N2 diagram
import argparse
parser = argparse.ArgumentParser()
parser.add_argument(
"--hide_visuals",
default=False,
action="store_true",
help="Do not show matplotlib figure or open N2 diagram in browser",
)
hide_viz = parser.parse_args().hide_visuals
# Setup the problem and run the analysis
prob = setup_problem(model=MissionAnalysisWithFuelBurn)
prob.run_model()
# Generate N2 diagram
om.n2(prob, outfile="minimal_integrator_n2.html", show_browser=not hide_viz)
# Create plot with results
import matplotlib.pyplot as plt
fig, axs = plt.subplots(2, 3, figsize=[9, 4.8], constrained_layout=True)
axs = axs.flatten() # change 2x3 mtx of axes into 4-element vector
# Define variables to plot
plot_vars = [
{"var": "fltcond|h", "name": "Altitude", "units": "ft"},
{"var": "fltcond|vs", "name": "Vertical speed", "units": "ft/min"},
{"var": "fltcond|Utrue", "name": "True airspeed", "units": "kn"},
{"var": "throttle", "name": "Throttle", "units": None},
{"var": "fuel_flow_calc.fuel_flow", "name": "Fuel flow", "units": "g/s"},
{"var": "weight", "name": "Weight", "units": "kg"},
]
for idx_fig, var in enumerate(plot_vars):
axs[idx_fig].set_xlabel("Range (nmi)")
axs[idx_fig].set_ylabel(f"{var['name']}" if var["units"] is None else f"{var['name']} ({var['units']})")
# Loop through each flight phase and plot the current variable from each
for phase in ["climb", "cruise", "descent"]:
axs[idx_fig].plot(
prob.get_val(f"mission.{phase}.range", units="nmi"),
prob.get_val(f"mission.{phase}.{var['var']}", units=var["units"]),
"-o",
c="tab:blue",
markersize=2.0,
)
fig.savefig("minimal_integrator_results.svg", transparent=True)
if not hide_viz:
plt.show()
The model should converge in a few iterations. The plot it generates should look like this:
The N2 diagram for the model is the following:
You can see that the weight is no longer constant. This results in a varying throttle in the cruise phase, unlike the constant throttle from the minimal example. Also notice that the fuel flow and throttle have the exact same shape, which makes sense because they are directly related by a factor of TSFC.
Summary
In this tutorial, we extended the previous minimal aircraft example to use an integrator to compute fuel burn. Our aircraft model is now an OpenMDAO group with a few more components in it. We compute fuel flow using thrust output by the aircraft model and TSFC. The integrator integrates fuel flow to compute fuel burn. Finally, we compute the aircraft weight by subtracting the fuel burned from the takeoff weight.
The time input for the integrator is connected automatically and the final integrated value from one phase is connected to the initial value for the following one with some Ben Brelje magic.
You’re encouraged to figure out how this works for yourself by looking at the source code for the PhaseGroup and TrajectoryGroup (these are used by the flight phase and mission analysis groups).
The final script looks like this:
"""
This example builds off the original minimal example,
but adds a numerical integrator to integrate fuel burn
and update the weight accordingly.
"""
# rst Imports (beg)
import openmdao.api as om
from openconcept.examples.minimal import Aircraft, setup_problem # build off this aircraft model
from openconcept.mission import BasicMission
from openconcept.utilities import Integrator
# rst Imports (end)
# rst Aircraft (beg)
class AircraftWithFuelBurn(om.Group):
"""
This model takes the simplified aircraft model from the minimal example, but adds
a fuel flow computation using TSFC and integrates it to compute fuel burn.
"""
# rst Options
def initialize(self):
self.options.declare("num_nodes", default=1, desc="Number of analysis points per phase")
self.options.declare("flight_phase", default=None) # required by OpenConcept but unused in this example
# rst Setup
def setup(self):
nn = self.options["num_nodes"]
# Add the aircraft model from the minimal example to build off of.
# Don't promote the weight from this model because we're going to compute a new
# one using the fuel burn.
# rst Simple aircraft (beg)
self.add_subsystem(
"simple_aircraft",
Aircraft(num_nodes=nn),
promotes_inputs=[
"fltcond|CL",
"throttle",
"fltcond|q",
"ac|geom|wing|S_ref",
"ac|weights|TOW",
"ac|propulsion|max_thrust",
"ac|aero|L_over_D",
],
promotes_outputs=["thrust", "drag"],
)
# rst Simple aircraft (end)
# Use an OpenMDAO ExecComp to compute the fuel flow rate using the thrust and TSFC
# rst Fuel flow (beg)
self.add_subsystem(
"fuel_flow_calc",
om.ExecComp(
"fuel_flow = TSFC * thrust",
fuel_flow={"units": "kg/s", "shape": nn},
TSFC={"units": "kg/N/s", "shape": 1},
thrust={"units": "N", "shape": nn},
has_diag_partials=True,
),
promotes_inputs=[("TSFC", "ac|propulsion|TSFC"), "thrust"],
)
# rst Fuel flow (end)
# Integrate the fuel flow rate to compute fuel burn
# rst Integrator (beg)
integ = self.add_subsystem(
"fuel_integrator", Integrator(num_nodes=nn, diff_units="s", time_setup="duration", method="simpson")
)
integ.add_integrand("fuel_burned", rate_name="fuel_flow", units="kg")
self.connect("fuel_flow_calc.fuel_flow", "fuel_integrator.fuel_flow")
# rst Integrator (end)
# Compute the current weight by subtracting the fuel burned from the takeoff weight
# rst Weight (beg)
self.add_subsystem(
"weight_calc",
om.ExecComp(
"weight = TOW - fuel_burned",
units="kg",
weight={"shape": nn},
TOW={"shape": 1},
fuel_burned={"shape": nn},
has_diag_partials=True,
),
promotes_inputs=[("TOW", "ac|weights|TOW")],
promotes_outputs=["weight"],
)
self.connect("fuel_integrator.fuel_burned", "weight_calc.fuel_burned")
# rst Weight (end)
# rst Mission (beg)
class MissionAnalysisWithFuelBurn(om.Group):
"""
OpenMDAO group for basic three-phase climb, cruise, descent mission.
The only top-level aircraft variable that the aircraft model uses is
the wing area, so that must be defined.
"""
def initialize(self):
self.options.declare("num_nodes", default=1, desc="Number of analysis points per phase")
def setup(self):
iv = self.add_subsystem("ac_vars", om.IndepVarComp(), promotes_outputs=["*"])
iv.add_output("ac|geom|wing|S_ref", val=25.0, units="m**2")
iv.add_output("ac|weights|TOW", val=5e3, units="kg")
iv.add_output("ac|propulsion|max_thrust", val=1e4, units="N")
iv.add_output("ac|propulsion|TSFC", val=20.0, units="g/kN/s")
iv.add_output("ac|aero|L_over_D", val=10.0)
# Define the mission
self.add_subsystem(
"mission",
BasicMission(aircraft_model=AircraftWithFuelBurn, num_nodes=self.options["num_nodes"]),
promotes_inputs=["ac|*"],
)
# rst Mission (end)
# rst Run (beg)
if __name__ == "__main__":
# Process command line argument to optionally not show figures and N2 diagram
import argparse
parser = argparse.ArgumentParser()
parser.add_argument(
"--hide_visuals",
default=False,
action="store_true",
help="Do not show matplotlib figure or open N2 diagram in browser",
)
hide_viz = parser.parse_args().hide_visuals
# Setup the problem and run the analysis
prob = setup_problem(model=MissionAnalysisWithFuelBurn)
prob.run_model()
# Generate N2 diagram
om.n2(prob, outfile="minimal_integrator_n2.html", show_browser=not hide_viz)
# Create plot with results
import matplotlib.pyplot as plt
fig, axs = plt.subplots(2, 3, figsize=[9, 4.8], constrained_layout=True)
axs = axs.flatten() # change 2x3 mtx of axes into 4-element vector
# Define variables to plot
plot_vars = [
{"var": "fltcond|h", "name": "Altitude", "units": "ft"},
{"var": "fltcond|vs", "name": "Vertical speed", "units": "ft/min"},
{"var": "fltcond|Utrue", "name": "True airspeed", "units": "kn"},
{"var": "throttle", "name": "Throttle", "units": None},
{"var": "fuel_flow_calc.fuel_flow", "name": "Fuel flow", "units": "g/s"},
{"var": "weight", "name": "Weight", "units": "kg"},
]
for idx_fig, var in enumerate(plot_vars):
axs[idx_fig].set_xlabel("Range (nmi)")
axs[idx_fig].set_ylabel(f"{var['name']}" if var["units"] is None else f"{var['name']} ({var['units']})")
# Loop through each flight phase and plot the current variable from each
for phase in ["climb", "cruise", "descent"]:
axs[idx_fig].plot(
prob.get_val(f"mission.{phase}.range", units="nmi"),
prob.get_val(f"mission.{phase}.{var['var']}", units=var["units"]),
"-o",
c="tab:blue",
markersize=2.0,
)
fig.savefig("minimal_integrator_results.svg", transparent=True)
if not hide_viz:
plt.show()
# rst Run (end)