Kepler’s Equation

Kepler’s Equation#

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

(3)#\[\begin{align} E - e \sin{E} = M \end{align}\]

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.

Using a BalanceComp and NewtonSolver#

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:

(4)#\[\begin{align} lhs(var) \cdot mult(var) = rhs(var) \end{align}\]

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.

BalanceComp also provides a way to supply the starting value for the implicit state variable (\(E\) in this case), via the guess_func argument. The supplied function should have a similar signature to the guess_nonlinear function of ImplicitComponent. 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 with a NewtonSolver is as follows:

Define a Group to contain the implicit system.
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

import openmdao.api as om

prob = om.Problem()

bal = om.BalanceComp()

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

# Use M (mean anomaly) as the initial guess for E (eccentric anomaly)
def guess_function(inputs, outputs, residuals):
    if np.abs(residuals['E']) > 1.0E-2:
        outputs['E'] = inputs['M']

bal.options['guess_func'] = guess_function

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

prob.model.set_input_defaults('M', 85.0, units='deg')

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

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

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

# Set up solvers
prob.model.linear_solver = om.DirectSolver()
prob.model.nonlinear_solver = om.NewtonSolver(solve_subsystems=False, maxiter=100, iprint=2)

prob.setup()

prob.set_val('M', 85.0, units='deg')
prob.set_val('E', 85.0, units='deg')
prob.set_val('ecc', 0.6)

prob.run_model()

M = prob.get_val('M')
E = prob.get_val('E')
print(f'M = {M}')
print(f'E = {E}')

NL: Newton 0 ; 0.3321557 1
NL: Newton 1 ; 0.117134648 0.352649822
NL: Newton 2 ; 0.00208780933 0.00628563452
NL: Newton 3 ; 7.36738647e-07 2.2180521e-06
NL: Newton 4 ; 9.19264664e-14 2.76757155e-13
NL: Newton Converged
M = [1.48352986]
E = [2.02317564]

Note

BalanceComp should be placed after those components which compute its inputs. Otherwise errors can arise due to timing issues.

Using Fixed-Point Iteration with an ImplicitComponent#

One approach to solving Kepler’s Equation is to use fixed-point iteration:

(5)#\[\begin{align} E_{n+1} = M + e \sin{E_{n}} \end{align}\]

We can achieve this behavior by defining a simple ImplicitComponent that “solves itself” by overriding its solve_nonlinear method.

class KeplerComp(om.ImplicitComponent):
    
    def initialize(self):
        self.options.declare('tol', types=float, default=1.0E-12, desc='convergence tolerance')
        self.options.declare('max_its', types=int, default=50, desc='maximum number of iterations')
    
    def setup(self):
        self.add_input('M', val=1.0, units='rad', desc='Mean anomaly')
        self.add_input('ecc', val=0.0, units=None, desc='eccentricity')
        self.add_output('E', val=1.0, units='rad', desc='Eccentric anomaly')
        
        self.declare_partials(of='E', wrt='*', method='cs')
        
    def apply_nonlinear(self, inputs, outputs, residuals):
        """ Compute the residual of output E. """
        M = inputs['M']
        e = inputs['ecc']
        E = outputs['E']
        residuals['E'] = E - M - e * np.sin(E)
        
    def solve_nonlinear(self, inputs, outputs):
        """ Determine the values of the implicit outputs that eliminate the residual. """
        M = inputs['M']
        e = inputs['ecc']
        E = outputs['E']
        
        # Compute the initial guess for iteration if the initial redisual is large.
        E = inputs['M'] if np.abs(E - M - e * np.sin(E)) > 1.0E-2 else outputs['E']
            
        for i in range(self.options['max_its']):
            E_old = E
            E = M + e * np.sin(E)
            if i == 0:
                print('iter ' + '   ' + 'max abs. error')
                print(5*'-' + '   '+ 14*'-')
            print(f'{i:>5} {np.max(np.abs(E - E_old)):16.12f}')
            if np.abs(E - E_old) < self.options['tol']:
                outputs['E'] = E
                break
        else:
            raise om.AnalysisError('Iteration Limit Exceeded')

Here the KeplerComp uses the solve_nonlinear method to drive it residual to zero through fixed point iteration.

Important

In OpenMDAO it doesn’t matter how the residuals of a system of eliminated. Inside solve_nonlinear you’re free to use a variety of techniques. As long as the resulting residual is reasonably close to zero, the derivatives across the implicit system will be accurate.

import numpy as np

import openmdao.api as om

prob = om.Problem()

prob.model.add_subsystem(name='kep_comp', subsys=KeplerComp(),
                         promotes_inputs=['M', 'ecc'],
                         promotes_outputs=['E'])


prob.setup()

prob.set_val('M', 85.0, units='deg')
prob.set_val('E', 1.0, units='deg')
prob.set_val('ecc', 0.6)

prob.run_model()
M = prob.get_val('M')
E = prob.get_val('E')
print(f'M = {M}')
print(f'E = {E}')

iter    max abs. error
-----   --------------
 0.597716818855
 0.074202079572
 0.020291240188
 0.005255936885
 0.001382786018
 0.000362353275
 0.000095052939
 0.000024927544
 0.000006537697
 0.000001714596
 0.000000449677
 0.000000117934
 0.000000030930
 0.000000008112
 0.000000002127
 0.000000000558
 0.000000000146
 0.000000000038
 0.000000000010
 0.000000000003
 0.000000000001
M = [1.48352986]
E = [2.02317564]

In contrast to the NewtonSolver approach, this approach roughly needs the following steps:

Define an ImplicitComponent where the apply_nonlinear method computes the residuals for Kepler’s Equation and solve_nonlinear find the value of the implicit output that eliminates the residuals.
Add the ImplicitComponent to the problem model.
Setup the problem, set values for the inputs, and run the model.

Which approach should you use?#

It’s really up to you. In this case, the Newton approach involves slightly more complicated System consisting of two components (the ExecComp and the BalanceComp) but it converges in just 4 iteratoins.

The ImplicitComponent involves only a single system, but writing the solve_nonlinear method with some error handling, initial guess generation, and print capability made for more coding on the part of the user.