dymos: A Python package for optimal control of multidisciplinary systems

Dymos is a library for optimizing control schedules for dynamic systems — sometimes referred to as optimal control or trajectory optimization. There are a number of other optimal control libraries that tackle similar kinds of problems, such as OTIS4 (Paris et al., 2006), GPOPSII (Patterson & Rao, 2014),and CASADI (Andersson et al., 2019). These tools all rely on gradient-based optimization to solve optimal control problems, though their methods of computing the gradients vary. Dymos is built on top of the OpenMDAO framework (Gray et al., 2019) and supports its modular derivative system which allows users to mix-and-match from finite-differencing, complex-step, hand-differentiated, and algorithmic differentiation. This flexibility allows Dymos to efficiently solve optimal control problems constructed with both ordinary differential equations (ODE) and differential-algebraic equations (DAE).


Summary
Dymos is a library for optimizing control schedules for dynamic systems -sometimes referred to as optimal control or trajectory optimization. There are a number of other optimal control libraries that tackle similar kinds of problems, such as OTIS4 (Paris et al., 2006), GPOPS-II (Patterson & Rao, 2014),and CASADI (Andersson et al., 2019). These tools all rely on gradient-based optimization to solve optimal control problems, though their methods of computing the gradients vary. Dymos is built on top of the OpenMDAO framework (Gray et al., 2019) and supports its modular derivative system which allows users to mix-and-match from finite-differencing, complex-step, hand-differentiated, and algorithmic differentiation. This flexibility allows Dymos to efficiently solve optimal control problems constructed with both ordinary differential equations (ODE) and differential-algebraic equations (DAE).
Dymos can also help solve more general optimization problems where dynamics are only one part in a larger system-level model with additional -potentially computationally expensivecalculations that come before and after the dynamic calculations. These broader problems are commonly referred to as co-design, controls-co-design, and multidisciplinary design optimization. Dymos provides specific APIs and features that make it possible to integrate traditional optimal-control models into a co-design context, while still supporting analytic derivatives that are necessary for computational efficiency in these complex use cases. An example of a co-design problem that was solved with Dymos is the coupled trajectory-thermal design of an electric vertical takeoff and landing aircraft where the thermal management and propulsion systems were designed simultaneously with the flight trajectories to ensure no components overheated (Hendricks et al., 2020).
In addition, there are constraints that typically need to be enforced: Path Constraints :p lb ≤ p (x, t,ū,d) ≤p ub In the mathematical sense what distinguishes optimal control from co-design is the particulars of which design variables and constraints are actually considered by the optimization. Pure optimal-control problems deal with a system of fixed design and seek to maximize performance by adjusting dynamic quantities (t,x,ū) such as position, speed, fuel-burned, and battery state-of-charge. Co-design problems simultaneously vary the static design parameters of a system (d) and its dynamic behavior (t,x,ū) to reach maximum performance.
In practice, the mathematical distinction is too rigid and a more practical distinction is made based on where the static and dynamic calculations are implemented and how complex each of them is. For very simple physical design parameters (e.g. the radius of a cannon ball, spring constants, linkage lengths, etc) it is common to integrate the design calculations directly into the ODE. Even though the calculations are static in nature, they can easily be coded as part of the ODE and still fit well into the optimal-control paradigm. The optimization structure thus looks like this: However, not all problems can be handled with such a compact implementation. For example if the physical design problem included shaping of an airfoil using expensive numerical solutions of partial differential equations (PDE) to predict drag, then one would not want to embed that PDE solver into the dynamic model. Instead the user could set up a coupled model with the PDE solver going first, and passing a table of data to be interpolated to the dynamic model. This effectively splits calculations up into static and dynamic components. This implementation structure is called co-design.
Traditionally, this co-design implementation would be done via sequential optimization with a manual outer design iteration between the static and dynamic models, potentially with different teams of people working on each one. One team would come up with a physical design using their own internal optimization setup. A second team takes the design and generates optimal-control profiles for it. Of course, the iterations do not need to be manual. It is also possible to set up an iterative loop around static and dynamic models to converge the problem numerically. A sequential co-design implementation looks like this: Dymos can support sequential co-design, but its unique value is that it also enables a more tightly-coupled co-design process with a single top level optimizer handling both parts of the problem simultaneously. Compared to sequential co-design, coupled co-design offers the potential to find better designs with much lower computational cost. However, it is also more challenging to implement because the top-level optimizer requires derivatives to be propagated between the static and dynamic parts of the model. Dymos overcomes this difficulty by providing APIs to exploit OpenMDAO's analytic derivative functionality at the model level. Data can be passed from the static model to the dynamic model and vice versa, allowing the construction of the coupled model for optimization.

ODE versus DAE
Optimal-control software typically requires that the dynamics of the system be defined as a set of ordinary differential equations (ODE) that use explicit functions to compute the rates of the state variables to be time-integrated. Sometimes the dynamics are instead posed as a set of differential-algebraic equations (DAE), where some residual equations need to be satisfied implicitly in order to solve for the state rates. From the perspective of an optimalcontrol or co-design problem both ODE and DAE formulations provide state rates that need to be integrated over time. The difference is that ODEs are explicit functions which are relatively easy to differentiate, but DAEs are implicit functions which are much more difficult to differentiate. Since the derivatives are needed to perform optimization, DAEs are more challenging to optimize.
One relatively common use case for DAEs is differential inclusions, in which the state trajectory is posed as a dynamic control and the traditional control schedule needed to achieve that trajectory is found using a nonlinear solver within the dynamic model (Seywald, 1994). For some problems this method provides a more natural and numerically-beneficial design space for the optimizer to traverse, but the nonlinear solver poses numerical challenges for computing derivatives for the optimizer. A simple approach to this is to just use finite-differences across the nonlinear solver, but this has been shown to be expensive and numerically unstable (Gray et al., 2014). Another option, taken by some optimal control libraries, is to apply monolithic algorithmic differentiation (Griewank, 2003) across the nonlinear solver. While it does provide accurate derivatives, the monolithic approach is expensive and uses a lot of memory (Kenway et al., 2019;Mader et al., 2008). The most efficient approach is to use a pair of analytic derivative approaches called the direct and adjoint methods, which were generalized in a single unified derivative equation (UDE) by Hwang and Martins (Hwang & Martins, 2018).
Dymos adopts the UDE approach, which uses a linear solver to compute total derivatives needed by the optimizer using only partial derivatives of the residual equations in the DAE. This approach offers two key advantages. First, partial derivatives of the DAE residual equations are much less computationally challenging to compute. Second, by using the OpenMDAO underpinnings of Dymos, users can construct their DAE in a modular fashion and combine various methods of computing the partial derivatives via finite-difference, complex-step (Martins et al., 2003), algorithmic differentiation, or hand differentiation as needed.

The Dymos perspective on optimal control
Dymos breaks the trajectory into portions of time called phases. Breaking the trajectory into phases provides several capabilities. Intermediate constraints along a trajectory can be enforced by applying a boundary constraint to a phase that begins or ends at the time of interest. For instance, the optimal trajectory of a launch vehicle may be required to ascend vertically to clear a launch tower before it pitches over on its way to orbit. Path constraints can be applied within each phase to bound some performance parameter within that phase. For example, reentry vehicles may need to adjust their trajectory to limit aerodynamic heating.
Each phase in a trajectory can use its own separate ODE. For instance, an aircraft with vertical takeoff and landing capability may use different ODEs for vertical flight and horizontal flight. ODEs are implemented as standard OpenMDAO models which are passed to phases at instantiation time with some additional annotations to identify the states, state-rates, and control inputs.
Every phase uses its own specific time discretization tailored to the dynamics in that portion of the trajectory. If one part of a trajectory has fast dynamics and another has slow dynamics, the time evolution can be broken into two phases with separate time discretizations.
In the optimal-control community there are a number of different techniques for discretizing the continuous optimal control problem into a form that can be solved by a nonlinear optimization algorithm; each one is called a transcription. Dymos supports two different collocation transcriptions: high-order Gauss-Lobatto (Herman & Conway, 1996) and Radau (Garg et al., 2009). Both of these represent state and control trajectories as piece-wise polynomials of at least 3rd order and are formulated in a way that makes it possible to efficiently compute the needed quantities to perform integration in a numerically rigorous fashion.
Dymos also allows the user to choose whether the optimization problem is solved using an explicit or implicit approach. Some caution with terminology must be taken here because the term "implicit" is often used to describe time integration schemes (e.g. backwards Euler), but that is not what is meant in an optimal-control context. Here, explicit propagation is one where the full state trajectory is computed starting from the initial value and propagating forward or from the final value and propagating backward. From the optimizer's perspective it will set values for the initial or final state (x), the design parameters (d), and the controls (ū) and can expect to be given a physically valid time evolution of the states as the output.
Wrapping an optimizer around an explicit propagation gives what is traditionally called a "shooting method" in the optimal-control world. In contrast, implicit propagation used within an optimization does not provide valid trajectories on its own. Instead, implicit methods add a discretized time-evolution of the state vector (x) as an additional design variable to the optimizer and add an associated set of defect constraints that must be driven to zero to enforce physics at some set of discrete points in time where the ODE is evaluated. The net effect is that the full state trajectory is only known once the optimization is fully converged.
In the context of the multidisciplinary design optimization field, explicit phases are similar to the multidisciplinary design feasible (MDF) optimization architecture and implicit phases are similar to the simultaneous analysis and design (SAND) optimization architecture (Martins & Lambe, 2013).
Both implicit and explicit phases are useful in different circumstances. Explicit propagation can seem to many like a more natural way to formulate the problem because it matches the way one would use time integration without optimization. However, when used with optimization explicit propagation is more computationally expensive, sensitive to singularities in the ODE, and potentially unable to converge to a valid solution. Implicit propagation tends to be less intuitive computationally, since it does not provide valid state histories without a converged optimization. The advantages of implicit propagation are that it tends to be faster, more numerically stable, and more scalable -though also highly sensitive to initial conditions and optimization scaling.
Dymos supports both explicit and implicit propagation for both its transcriptions, and even allows mixtures of implicitly and explicitly propagated states within a phase. This flexibility is valuable because it allows users to tailor their optimization to suit their needs. Switching transcriptions and changing from implicit to explicit requires very minor code changestypically a single line in the run-script. Examples of how to swap between them are given in the code sample below.

Choice of optimization algorithm
Dymos is not distributed with an optimizer, but relies on the optimizers that are available in the OpenMDAO installation. OpenMDAO ships with an interface to the optimizers in SciPy (Virtanen et al., 2020), and an additional wrapper for the pyoptsparse library (Wu et al., 2020) which has more powerful optimizer options such as SNOPT (Gill et al., 2005) and IPOPT (Wächter & Biegler, 2006). OpenMDAO also allows users to integrate their own optimizer of choice, which Dymos can then seamlessly use with without any additional modifications. For simple problems, Scipy's SLSQP optimizer generally works fine. On more challenging optimal-control problems higher-quality optimizers are important for getting good performance.
Though one could technically choose any optimization algorithm, Dymos is designed to work primarily with gradient-based algorithms. In general, optimal-control and co-design problems will have both a very large number of design variables and a very large number of constraints. Both of these issues make gradient-based methods the strongly-preferred choice. Gradientfree methods could potentially be used in certain narrow circumstances with problems built using purely explicit phases and set up intentionally to have a small set of design variables and constraints.

Statement of Need
When dealing with the design of complex systems that include transient behavior, co-design becomes critical (Garcia-Sanz, 2019). Broadly there are two approaches: sequential co-design or coupled co-design (Fathy et al., 2001;Peters et al., 2009). The best choice depends on the degree of interaction, or coupling, between various sub-systems. If the coupling is strong a coupled co-design approach is necessary to achieve the best performance.
Though there are a number of effective optimal-control libraries, they tend to assume that they are on top of the modeling stack. They frame every optimization problem as if it were a pure optimal-control problem, and hence are best suited to be used in a sequential co-design style. This poses large challenges when expanding to tightly-coupled problems, where the interactions between the static and dynamic systems are very strong.
Dymos provides a set of unique capabilities that make coupled co-design possible via efficient gradient-based optimization methods. It provides differentiated time-integration schemes that can generate transient models from user provided ODEs, along with APIs that enable users to couple these transient models with other models to form the co-design system while carrying the differentiation through that coupling. It also supports efficient differentiation of DAEs that include implicit relationships, which allows for a much broader set of possible ways to pose transient models. These two features combined make Dymos capable of handling coupled co-design problems in a manner that is more efficient than a pure optimal-control approach.

Selected applications of Dymos
Dymos has been used to demonstrate the coupling of flight dynamics and subsystem thermal constraints in electrical aircraft applications (Falck et al., 2017;Hendricks et al., 2020). NASA's X-57 "Maxwell" is using Dymos for mission planning to maximize data collection while abiding the limits of battery storage capacity and subsystem temperatures (Schnulo et al., 2018(Schnulo et al., , 2019. Other authors have used Dymos to perform studies of aircraft acoustics  and the design of supersonic aircraft with thermal fuel management systems (Jasa et al., 2018).

Optimal-control example: Brachistochrone
As a simple use case of Dymos, consider the classic brachistochrone optimal-control problem. There is a bead sliding along a frictionless wire strung between two points of different heights, and we seek the shape of the wire such that the bead travels from start to finish in the shortest time. The time-varying control is the angle of the wire at each point in time and there are no design parameters, which makes this a pure optimal-control problem. The built-in plotting utility in Dymos will plot all relevant quantities vs time: The more traditional way to view the brachistochrone solution is to view the actual shape of the wire (i.e. y vs x):

Coupled co-design example: Designing a cannonball
This co-design example seeks to find the best size cannonball to maximize range, considering aerodynamic drag subject to a limit on initial kinetic energy. Given the same kinetic energy, a lighter ball will go faster, and hence farther, if aerodynamic drag is ignored. Heavier cannonballs will have more inertia to counteract drag. There is a balance between these two effects, which the optimizer seeks to find.
Here the static calculations are to find the mass and frontal area of the cannonball, given its radius. Then the ODE takes the mass and area as inputs and via Dymos can compute the total range. For demonstration purposes the trajectory is broken up into an ascent and descent phase, with the break being set up exactly at the apogee of the flight path.