# 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:

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:

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.