Mission Solver Code Structure
This is a high level overview of how the mission solver functions. The purpose is to show the structure that is used for an existing mission, and show where changes should be made if different functionality is desired.
File Structure
Mission scripts are split into two folders in the SUAVE repository. The first is in trunk/SUAVE/Analyses/Mission/Segments, and the second is in trunk/SUAVE/Methods/Missions/Segments. As with other types of analyses and methods, the distinction between these is that the Analyses folder contains classes that are built to use functions stored in the Methods folder. This division is done to make it easier to build new analysis classes using a mix of available methods.
A typical mission segment analysis file contains four keys parts. The first specifies default user inputs, unknowns, and residuals. The inputs are used to provide the analysis with conditions that need to be met, while the unknowns and residuals are used as part of the solution process. The second sets the initialization functions for the analysis, which are run at the beginning. The third picks the convergence method and specifies the functions that will be used during iteration. The fourth finalizes the data and processes it for results output.
Initialization
For this tutorial, we will be considering the constant speed constant altitude cruise segment. The files are available here (Analysis) and here (Method). This class also inherits information from more general segment classes, which include many of the processing functions. As with other segments, the user will specify key conditions. For this case, altitude, air speed, and distance are the necessary inputs. If the user does not specify an altitude, it will be taken automatically from the last value in the previous segment. These inputs must be specified in some way for the mission segment to be evaluated. They are shown below as well:
self.altitude = None
self.air_speed = 10. * Units['km/hr']
self.distance = 10. * Units.km
The other set of segment specific initial values are the values used for solving the segment (typically this means satisfying a force balance at every evaluation point). These can be changed by the user if needed, but the default values should perform fine for most cases.
self.state.unknowns.throttle = ones_row(1) * 0.5
self.state.unknowns.body_angle = ones_row(1) * 0.0
self.state.residuals.forces = ones_row(2) * 0.0
Here throttle and body angle are the unknowns, and the values shown here are the values they will start at. The residuals will be computed based on these unknowns, so their initial value is not important. Instead they are initialized just to create the necessary data structure. The ones_row line will create a numpy array with the number of elements needed for evaluation.
Evaluation Details
Most of the missions in SUAVE, including this one, are broken into several points in time based on a Chebyshev polynomial. This causes the points to be closer together at either end of the segment. The choice of a Chebyshev polynomial (which creates cosine spacing) provides better convergence and smoothness properties versus other methods such as linear spacing.
At each of these points the aerodynamic analysis is queried to find CL and CD, which are then converted to lift and drag. These values will be dependent on the body angle unknown and other aerodynamic parameters. Thrust is found from the vehicle’s energy network, which is dependent on the throttle unknown. A weight is determined by looking at the initial weight and subsequent mass rate (typically corresponding with fuel burn). In this cruise segment, these forces are summed in 2D and the results are put in the residuals. The functions needed to arrive these forces are found in the Update Conditions section of the Analysis file. This section is also shown below in one of the steps to create a new mission.
Once the evaluation process has been performed at all points, the unknowns and residuals are fed back to the solve routine, which in this case is scipy’s fsolve. The file that performs this process is here. This routine continues evaluating the points until convergence is reached. Once this happens, post processing is done to put the data in the results output.
Using Multiple Segments
Multiple segments can be run sequentially by appending them in the desired order. Examples of this are in all the tutorial files that have an aircraft fly a full mission. In addition, the full mission can be run simultaneously will all segment constraints used together. If you are interested in doing something like this, please ask us about it on our forum.
Process Summary
Mission Setup
- Initializes default values for unknowns
- Initializes set of functions used to determine residuals
- Reads user input for segment parameters
- Adds the analysis group to be used (including the vehicle and items like atmosphere)
- Appends segments in order
Evaluate
- Varies unknowns until residual convergence is reached using scipy’s fsolve
- Repeats process for each segment until full mission is complete
Adding New Mission Segments
The segment described above uses two unknowns to solve force residuals in two dimensions. This general setup works well for many problems of interest, but SUAVE is designed to accommodate other mission analysis types as well. A user may want to add control surface deflection and solve for moments as well, or look at forces in all three dimensions.
In addition, a user may want to modify how the mission is flown, as is done with the many other segments currently available. They may want to modify how the mission is solved, such as is done in our single point evaluation segments where finite differencing is not relevant.
Here we will explain the process of modifying our constant speed constant rate climb segment to be constant throttle constant speed. This still uses 2D force balance but changes the profile. There are four functions that are modified here. The first is shown below. The functions can be found in here and here
def initialize_conditions(segment,state):
# unpack
climb_rate = segment.climb_rate
air_speed = segment.air_speed
alt0 = segment.altitude_start
altf = segment.altitude_end
t_nondim = state.numerics.dimensionless.control_points
conditions = state.conditions
# check for initial altitude
if alt0 is None:
if not state.initials: raise AttributeError('initial altitude not set')
alt0 = -1.0 * state.initials.conditions.frames.inertial.position_vector[-1,2]
# discretize on altitude
alt = t_nondim * (altf-alt0) + alt0
# process velocity vector
v_mag = air_speed
v_z = -climb_rate # z points down
v_x = np.sqrt( v_mag**2 - v_z**2 )
# pack conditions
conditions.frames.inertial.velocity_vector[:,0] = v_x
conditions.frames.inertial.velocity_vector[:,2] = v_z
conditions.frames.inertial.position_vector[:,2] = -alt[:,0] # z points down
conditions.freestream.altitude[:,0] = alt[:,0] # positive altitude in this context
This function initializes speed and altitude based on the given climb rate, airspeed, and altitude end points. t_nondim gives nondimensional time in cosine spacing from 0 to 1 in order to pick the values at the points to be evaluated. Unfortunately, when we use constant throttle we cannot know beforehand exactly how altitude (or climb rate in this case) will vary with time, so altitude cannot be spaced with this method. Instead a different function is used to initialize conditions:
def initialize_conditions(segment,state):
# unpack
throttle = segment.throttle
air_speed = segment.air_speed
alt0 = segment.altitude_start
altf = segment.altitude_end
t_nondim = state.numerics.dimensionless.control_points
conditions = state.conditions
# check for initial altitude
if alt0 is None:
if not state.initials: raise AttributeError('initial altitude not set')
alt0 = -1.0 * state.initials.conditions.frames.inertial.position_vector[-1,2]
# pack conditions
conditions.propulsion.throttle[:,0] = throttle
conditions.frames.inertial.velocity_vector[:,0] = air_speed # start up value
Here only the throttle and air speed are loaded in, and discretization of other values will need to occur later so that it is part of the iteration loop. This requires a new function that updates the altitude differentials.
def update_differentials_altitude(segment,state):
# unpack
t = state.numerics.dimensionless.control_points
D = state.numerics.dimensionless.differentiate
I = state.numerics.dimensionless.integrate
# Unpack segment initials
alt0 = segment.altitude_start
altf = segment.altitude_end
conditions = state.conditions
r = state.conditions.frames.inertial.position_vector
v = state.conditions.frames.inertial.velocity_vector
# check for initial altitude
if alt0 is None:
if not state.initials: raise AttributeError('initial altitude not set')
alt0 = -1.0 * state.initials.conditions.frames.inertial.position_vector[-1,2]
# get overall time step
vz = -v[:,2,None] # Inertial velocity is z down
dz = altf- alt0
dt = dz / np.dot(I[-1,:],vz)[-1] # maintain column array
# Integrate vz to get altitudes
alt = alt0 + np.dot(I*dt,vz)
# rescale operators
t = t * dt
# pack
t_initial = state.conditions.frames.inertial.time[0,0]
state.conditions.frames.inertial.time[:,0] = t_initial + t[:,0]
conditions.frames.inertial.position_vector[:,2] = -alt[:,0] # z points down
conditions.freestream.altitude[:,0] = alt[:,0] # positive altitude in this context
return
In this function, t, D, and I are numpy arrays that allow approximate differentiation and integration. Since the total time is not known without determining the climb rate, we must first determine the time required to reach the final altitude. The line dt = dz / np.dot(I[-1,:],vz)[-1] does this with the integrator providing the amount of altitude gained if the velocities were spread across just one second instead of the full segment time. This gives the scaling quantity dt that is then used to get the altitude at every point in alt = alt0 + np.dot(I*dt,vz). The values for altitude are then are then packed for use in other functions.
The above allows us to deal with discretization without a known profile, but we also must calculate the velocity in order to use this. This is done with another added function.
def update_velocity_vector_from_wind_angle(segment,state):
# unpack
conditions = state.conditions
v_mag = segment.air_speed
alpha = state.unknowns.wind_angle[:,0][:,None]
theta = state.unknowns.body_angle[:,0][:,None]
# Flight path angle
gamma = theta-alpha
# process
v_x = v_mag * np.cos(gamma)
v_z = -v_mag * np.sin(gamma) # z points down
# pack
conditions.frames.inertial.velocity_vector[:,0] = v_x[:,0]
conditions.frames.inertial.velocity_vector[:,2] = v_z[:,0]
return conditions
This uses our new set of unknowns to determine the velocities.
Additionally, since the unknowns are different we must change the function that unpacks them. Wind angle does not need to be stored so it is not included here.
def unpack_body_angle(segment,state):
# unpack unknowns
theta = state.unknowns.body_angle
# apply unknowns
state.conditions.frames.body.inertial_rotations[:,1] = theta[:,0]
We now add these functions to the segment process list.
# --------------------------------------------------------------
# Initialize - before iteration
# --------------------------------------------------------------
initialize = self.process.initialize
initialize.expand_state = Methods.expand_state
initialize.differentials = Methods.Common.Numerics.initialize_differentials_dimensionless
initialize.conditions = Methods.Climb.Constant_Throttle_Constant_Speed.initialize_conditions
initialize.velocities = Methods.Climb.Constant_Throttle_Constant_Speed.update_velocity_vector_from_wind_angle
initialize.differentials_altitude = Methods.Climb.Constant_Throttle_Constant_Speed.update_differentials_altitude
and
# Unpack Unknowns
iterate.unknowns = Process()
iterate.unknowns.mission = Methods.Climb.Constant_Throttle_Constant_Speed.unpack_body_angle
# Update Conditions
iterate.conditions = Process()
iterate.conditions.velocities = Methods.Climb.Constant_Throttle_Constant_Speed.update_velocity_vector_from_wind_angle
iterate.conditions.differentials_a = Methods.Climb.Constant_Throttle_Constant_Speed.update_differentials_altitude
iterate.conditions.differentials_b = Methods.Common.Numerics.update_differentials_time
iterate.conditions.acceleration = Methods.Common.Frames.update_acceleration
iterate.conditions.altitude = Methods.Common.Aerodynamics.update_altitude
iterate.conditions.atmosphere = Methods.Common.Aerodynamics.update_atmosphere
iterate.conditions.gravity = Methods.Common.Weights.update_gravity
iterate.conditions.freestream = Methods.Common.Aerodynamics.update_freestream
iterate.conditions.orientations = Methods.Common.Frames.update_orientations
iterate.conditions.aerodynamics = Methods.Common.Aerodynamics.update_aerodynamics
iterate.conditions.stability = Methods.Common.Aerodynamics.update_stability
iterate.conditions.propulsion = Methods.Common.Energy.update_thrust
iterate.conditions.weights = Methods.Common.Weights.update_weights
iterate.conditions.forces = Methods.Common.Frames.update_forces
iterate.conditions.planet_position = Methods.Common.Frames.update_planet_position
If you have any questions that are not answered in other tutorials or the FAQ please ask on our forum page. This is also the place to go if you want help building a more elaborate evaluation, such as one that includes moments.