<jupyter_start><jupyter_text>Kepler’s EquationThis example will demonstrate the use of OpenMDAO for solving an implicit equation commonly found in astrodynamics, Kepler’s Equation:\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 NewtonSolverIn 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:\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:1. Define a Group to contain the implicit system.2. To that Group, add components which provide, $M$, $e$, and the left-hand side of Kepler’s equation.3. Add a linear and nonlinear solver to the Group, since the default solvers do not iterate.4. Setup the problem, set values for the inputs, and run the model.<jupyter_code>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.add_subsystem(name='balance', subsys=bal,
                         promotes_inputs=['M'],
                         promotes_outputs=['E'])

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

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')

# 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}')
# Check that results are correct.

from openmdao.utils.assert_utils import assert_near_equal

assert_near_equal(prob.get_val('E'), 2.02317564, tolerance=1.0E-6)<jupyter_output><empty_output><jupyter_text>Using Fixed-Point Iteration with an ImplicitComponentOne approach to solving Kepler's Equation is to use fixed-point iteration:\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.<jupyter_code>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')<jupyter_output><empty_output><jupyter_text>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.```<jupyter_code>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}')
# Check that results are correct.

from openmdao.utils.assert_utils import assert_near_equal

assert_near_equal(prob.get_val('E'), 2.02317564, tolerance=1.0E-6)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12297397769516728, 'fraction_numerical': 0.011895910780669145, 'mean_word_length': 3.481012658227848, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '260578', 'license': 'Apache-2.0', 'n_tokens_mistral': 2195, 'n_tokens_neox': 2068, 'n_words': 681, 'path': 'openmdao/docs/openmdao_book/examples/keplers_equation.ipynb', 'repo_name': 'sebasgo/OpenMDAO', 'size': '6725'}