Kepler’s Equation Example - Solving an Implicit Equaiton

This example will demonstrate the use of OpenMDAO for solving an implicit equation commonly found in astrodynamics, Kepler’s Equation:

\[E - e \sin{E} = M\]

Here \(M\) is the mean anomaly, \(E\) is the eccentric anomaly, and \(e\) is the eccentricity of the orbit.

If we know the eccentric anomaly, computing the mean anomaly is trivial. However, solving for the eccentric anomaly when given the mean anomaly must be done numerically. We’ll do so using a nonlinear solver. In OpenMDAO, solvers converge all implicit state variables in a Group by driving their residuals to zero.

In an effort to simplify things for users, OpenMDAO features a Balance component. For each implicit state variable we assign to the balance, it solves the following equation:

\[lhs(var) \cdot mult(var) = rhs(var)\]

The \(mult\) term is an optional multiplier than can be applied to the left-hand side (LHS) of the equation. For our example, we will assign the right-hand side (RHS) to the mean anomaly (\(M\)), and the left-hand side to \(E - e \sin{E}\)

In this implementation, we rely on an ExecComp to compute the value of the LHS.

Implicit components in OpenMDAO also provide an “initial guess” method, guess_nonlinear, which provides the starting value for the implicit state variable (\(E\) in this case) for the nonlinear solver. When solving Kepler’s equation, using \(M\) as the initial guess for \(E\) is a good starting point.

In summary, the recipe for solving Kepler’s equation is as follows:

  • Define a problem with a Group as its model.
  • To that Group, add components which provide, \(M\), \(e\), and the left-hand side of Kepler’s equation.
  • Add a linear and nonlinear solver to the Group, since the default solvers do not iterate.
  • Setup the problem, set values for the inputs, and run the model.
import numpy as np
from numpy.testing import assert_almost_equal

from openmdao.api import Problem, Group, IndepVarComp, BalanceComp,    ExecComp, DirectSolver, NewtonSolver

prob = Problem(model=Group())

ivc = IndepVarComp()

ivc.add_output(name='M',
               val=0.0,
               units='deg',
               desc='Mean anomaly')

ivc.add_output(name='ecc',
               val=0.0,
               units=None,
               desc='orbit eccentricity')

bal = BalanceComp()

bal.add_balance(name='E', val=0.0, units='rad', eq_units='rad', rhs_name='M')

# Override the guess_nonlinear method, always initialize E to the value of M
def guess_func(inputs, outputs, residuals):
    outputs['E'] = inputs['M']

bal.guess_nonlinear = guess_func

# ExecComp used to compute the LHS of Kepler's equation.
lhs_comp = ExecComp('lhs=E - ecc * sin(E)',
                    lhs={'value': 0.0, 'units': 'rad'},
                    E={'value': 0.0, 'units': 'rad'},
                    ecc={'value': 0.0})

prob.model.add_subsystem(name='ivc', subsys=ivc, promotes_outputs=['M', 'ecc'])

prob.model.add_subsystem(name='balance', subsys=bal,
                         promotes_inputs=['M'],
                         promotes_outputs=['E'])

prob.model.add_subsystem(name='lhs_comp', subsys=lhs_comp,
                         promotes_inputs=['E', 'ecc'])

# Explicit connections
prob.model.connect('lhs_comp.lhs', 'balance.lhs:E')

# Setup solvers
prob.model.linear_solver = DirectSolver()
prob.model.nonlinear_solver = NewtonSolver()
prob.model.nonlinear_solver.options['maxiter'] = 100
prob.model.nonlinear_solver.options['iprint'] = 0

prob.setup()

prob['M'] = 85.0
prob['ecc'] = 0.6

prob.run_model()

print(np.degrees(prob['E']))
[ 115.91942563]
print('E (deg) = ', np.degrees(prob['E'][0]))
E (deg) =  115.91942563

Summary

We built a model consisting of a Group in which the nonlinear solver solves Kepler’s equation through the use of a BalanceComp. We also demonstrated the use of the guess_nonlinear method in implicit components.