Implicit Relationships

It is possible to compute outputs implicitly by defining a residual variable in terms of the output and inputs.

In the first example, we solve a quadratic equation. This quadratic has two solutions: 1 and 3. Depending on the starting value of the output variable, OpenMDAO will find one root or the other.

from openmdao.api import Problem
from openmdao.api import ScipyKrylov, NewtonSolver, NonlinearBlockGS
from omtools.api import Group, ImplicitComponent
import numpy as np


class ExampleApplyNonlinear(ImplicitComponent):
    def setup(self):
        g = self.group

        with g.create_group('sys') as group:
            group.create_indep_var('a', val=1)
            group.create_indep_var('b', val=-4)
            group.create_indep_var('c', val=3)
        a = g.declare_input('a')
        b = g.declare_input('b')
        c = g.declare_input('c')

        x = g.create_implicit_output('x')
        y = a * x**2 + b * x + c

        x.define_residual(y)
        self.linear_solver = ScipyKrylov()
        self.nonlinear_solver = NewtonSolver(solve_subsystems=False)


prob = Problem()
prob.model = Group()
prob.model.add_subsystem('example', ExampleApplyNonlinear())
prob.setup(force_alloc_complex=True)
prob.run_model()

print('using default x=1')
import omtools.examples.valid.ex_implicit_apply_nonlinear as ex

print('setting x=1.9')
ex.prob.set_val('x', 1.9)
ex.prob.run_model()
print(ex.prob['x'])

print('setting x=2.1')
ex.prob.set_val('x', 2.1)
ex.prob.run_model()
print(ex.prob['x'])

using default x=1

============
comp_example
============
NL: Newton Converged in 0 iterations
setting x=1.9

============
comp_example
============
NL: Newton Converged in 7 iterations
[1.]
setting x=2.1

============
comp_example
============
NL: Newton Converged in 7 iterations
[3.]

The expressions for the residuals will tell OpenMDAO to construct the relevant Component objects, but they will be part of a Problem within the generated ImplicitComponent object. As a result, the expressions defined above do not translate to Component objects in the outer Problem whose model is displayed in the n2 diagram below.

For especially complicated problems, where the residual may converge for multiple solutions, or where the residual is difficult to converge over some interval, omtools provides an API for bracketing solutions.

from openmdao.api import Problem
from openmdao.api import ScipyKrylov, NewtonSolver, NonlinearBlockGS
from omtools.api import Group, ImplicitComponent
import numpy as np


class ExampleBracketedScalar(ImplicitComponent):
    def setup(self):
        g = self.group

        with g.create_group('sys') as group:
            group.create_indep_var('a', val=1)
            group.create_indep_var('b', val=-4)
            group.create_indep_var('c', val=3)
        a = g.declare_input('a')
        b = g.declare_input('b')
        c = g.declare_input('c')

        x = g.create_implicit_output('x')
        y = a * x**2 + b * x + c

        x.define_residual_bracketed(
            y,
            x1=0,
            x2=2,
        )


prob = Problem()
prob.model = Group()
prob.model.add_subsystem('example', ExampleBracketedScalar())
prob.setup(force_alloc_complex=True)
prob.run_model()

print('x', prob['x'].shape)
print(prob['x'])

x (1,)
[1.]

Brackets may also be specified for multidimensional array values.

from openmdao.api import Problem
from openmdao.api import ScipyKrylov, NewtonSolver, NonlinearBlockGS
from omtools.api import Group, ImplicitComponent
import numpy as np


class ExampleBracketedArray(ImplicitComponent):
    def setup(self):
        g = self.group

        with g.create_group('sys') as group:
            group.create_indep_var('a', val=[1, -1])
            group.create_indep_var('b', val=[-4, 4])
            group.create_indep_var('c', val=[3, -3])
        a = g.declare_input('a', shape=(2, ))
        b = g.declare_input('b', shape=(2, ))
        c = g.declare_input('c', shape=(2, ))

        x = g.create_implicit_output('x', shape=(2, ))
        y = a * x**2 + b * x + c

        x.define_residual_bracketed(
            y,
            x1=[0, 2.],
            x2=[2, np.pi],
        )


prob = Problem()
prob.model = Group()
prob.model.add_subsystem('example', ExampleBracketedArray())
prob.setup(force_alloc_complex=True)
prob.run_model()

print('x', prob['x'].shape)
print(prob['x'])

x (2,)
[1. 3.]