OpenMDAO implicit state with finite differences

Question

I would like to use an implicit variable and a nonlinear solver to iterate on its value in a system using finite differences for derivative calculation. I understand the Newton solver needs analytical derivatives but the non-linear Gauss-Siedel does not seem to be iterating correctly either. To be sure that I had coded things up appropriately, I tested this approach on the intersect_parabola_line.py example provided in the package.

from __future__ import print_function
from openmdao.api import Component, Group, Problem, Newton, ScipyGMRES,NLGaussSeidel

class Line(Component):
    """Evaluates y = -2x + 4."""
    def __init__(self):
        super(Line, self).__init__()
        self.add_param('x', 1.0)
        self.add_output('y', 0.0)
        # User can change these.
        self.slope = -2.0
        self.intercept = 4.0

    def solve_nonlinear(self, params, unknowns, resids):
        """ y = -2x + 4 """
        x = params['x']
        m = self.slope
        b = self.intercept
        unknowns['y'] = m*x + b

class Parabola(Component):
    """Evaluates y = 3x^2 - 5"""
    def __init__(self):
        super(Parabola, self).__init__()
        self.add_param('x', 1.0)
        self.add_output('y', 0.0)
        # User can change these.
        self.a = 3.0
        self.b = 0.0
        self.c = -5.0

    def solve_nonlinear(self, params, unknowns, resids):
        """ y = 3x^2 - 5 """
        x = params['x']
        a = self.a
        b = self.b
        c = self.c
        unknowns['y'] = a*x**2 + b*x + c

class Balance(Component):
    """Evaluates the residual y1-y2"""
    def __init__(self):
        super(Balance, self).__init__()
        self.add_param('y1', 0.0)
        self.add_param('y2', 0.0)
        self.add_state('x', 5.0)

    def solve_nonlinear(self, params, unknowns, resids):
        """This component does no calculation on its own. It mainly holds the
        initial value of the state. An OpenMDAO solver outside of this
        component varies it to drive the residual to zero."""
        pass

    def apply_nonlinear(self, params, unknowns, resids):
        """ Report the residual y1-y2 """
        y1 = params['y1']
        y2 = params['y2']
        resids['x'] = y1 - y2

if __name__ == '__main__':
    top = Problem()
    root = top.root = Group()
    root.add('line', Line())
    root.add('parabola', Parabola())
    root.add('bal', Balance())

    root.connect('line.y', 'bal.y1')
    root.connect('parabola.y', 'bal.y2')
    root.connect('bal.x', 'line.x')
    root.connect('bal.x', 'parabola.x')
    root.deriv_options['type'] = 'fd'

    root.nl_solver = NLGaussSeidel() #Newton()
    root.ln_solver = ScipyGMRES()
    root.nl_solver.options['iprint'] = 2

    top.setup()

    # Positive solution
    top['bal.x'] = 7.0
    root.list_states()
    top.run()
    print('Positive Solution x=%f, line.y=%f, parabola.y=%f' % (top['bal.x'], top['line.y'], top['parabola.y']))

    # Negative solution
    top['bal.x'] = -7.0
    root.list_states()
    top.run()
    print('Negative Solution x=%f, line.y=%f, parabola.y=%f' % (top['bal.x'], top['line.y'], top['parabola.y']))

The output that I get is:

States in model:

bal.x
Value: 7.0
Residual: 0.0

[root] NL: NLN_GS   1 | 152 1
[root] NL: NLN_GS   2 | 152 1
[root] NL: NLN_GS   2 | Converged in 2 iterations
Positive Solution x=7.000000, line.y=-10.000000, parabola.y=142.000000

States in model:

bal.x
Value: -7.0
Residual: -152.0

[root] NL: NLN_GS   1 | 124 1
[root] NL: NLN_GS   2 | 124 1
[root] NL: NLN_GS   2 | Converged in 2 iterations
Negative Solution x=-7.000000, line.y=18.000000, parabola.y=142.000000

Any tips would be appreciated. I am using Python 2.7.13 and OpenMDAO 1.7.3 on OSX.

Kenneth Moore Kenneth Moore · Accepted Answer · 2017-11-09T17:02:19

So, Nonlinear Gauss Seidel is not really going to be able to converge the model with a balance component and implicit state. You should really use Newton for that. To make that work, you want to also make sure that the models that Newton converges use finite difference too. When you set fd in the root group, that means that derivative calculations that originate above that system see an approximated fd on that whole system.

To make this work, try this:

top = Problem()
root = top.root = Group()
comp1 = root.add('line', Line())
comp2 = root.add('parabola', Parabola())
comp3 = root.add('bal', Balance())

root.connect('line.y', 'bal.y1')
root.connect('parabola.y', 'bal.y2')
root.connect('bal.x', 'line.x')
root.connect('bal.x', 'parabola.x')
root.deriv_options['type'] = 'fd'
comp1.deriv_options['type'] = 'fd'
comp2.deriv_options['type'] = 'fd'
comp3.deriv_options['type'] = 'fd'

OpenMDAO implicit state with finite differences

1 Answers