ImplicitComponent

ImplicitComponent#

Implicit variables are those that are computed as an implicit function of other variables. For instance, \(y\) would be an implicit variable, given that it is computed by solving \(cos(x⋅y)−z⋅y=0\). In OpenMDAO, implicit variables are defined as the outputs of components that inherit from ImplicitComponent.

In the above implicit expression, \(y\) is the implicit variable while \(x\) and \(z\) would be considered inputs.

ImplicitComponent Methods#

The implementation of each method will be illustrated using a simple implicit component that computes the output \(x\) implicitly via a quadratic equation, \(ax^2+bx+c=0\), where \(a\), \(b\), and \(c\) are inputs to the component.

import openmdao.api as om


class QuadraticComp(om.ImplicitComponent):
    """
    A Simple Implicit Component representing a Quadratic Equation.

    R(a, b, c, x) = ax^2 + bx + c

    Solution via Quadratic Formula:
    x = (-b + sqrt(b^2 - 4ac)) / 2a
    """

setup() :

Declare input and output variables via add_input and add_output. Information like variable names, sizes, units, and bounds are declared. Also, declare partial derivatives that this component provides. Here we use the wild card to say that this component provides derivatives of all implicit residuals with respect to all inputs and outputs.

def setup(self):
    self.add_input('a', val=1., tags=['tag_a'])
    self.add_input('b', val=1.)
    self.add_input('c', val=1.)
    self.add_output('x', val=0., tags=['tag_x'])

setup_partials():

Declare partial derivatives that this component provides. Here we use the wild card to say that this component provides derivatives of all implicit residuals with respect to all inputs and outputs.

def setup_partials(self):
    self.declare_partials(of='*', wrt='*')

apply_nonlinear(inputs, outputs, residuals) :

Compute the residuals, given the inputs and outputs.

def apply_nonlinear(self, inputs, outputs, residuals):
    a = inputs['a']
    b = inputs['b']
    c = inputs['c']
    x = outputs['x']
    residuals['x'] = a * x ** 2 + b * x + c

solve_nonlinear(inputs, outputs) :

[Optional] Compute the outputs, given the inputs.

def solve_nonlinear(self, inputs, outputs):
    a = inputs['a']
    b = inputs['b']
    c = inputs['c']
    outputs['x'] = (-b + (b ** 2 - 4 * a * c) ** 0.5) / (2 * a)

linearize(inputs, outputs, partials) :

[Optional] The component’s partial derivatives of interest are those of the residuals with respect to the inputs and the outputs. If the user computes the partial derivatives explicitly, they are provided here. If the user wants to implement partial derivatives in a matrix-free way, this method provides a place to perform any necessary assembly or pre-processing for the matrix-vector products. Regardless of how the partial derivatives are computed, this method provides a place to perform any relevant factorizations for directly solving or preconditioning the linear system.

def linearize(self, inputs, outputs, partials):
    a = inputs['a']
    b = inputs['b']
    c = inputs['c']
    x = outputs['x']

    partials['x', 'a'] = x ** 2
    partials['x', 'b'] = x
    partials['x', 'c'] = 1.0
    partials['x', 'x'] = 2 * a * x + b

    self.inv_jac = 1.0 / (2 * a * x + b)

apply_linear(inputs, outputs, d_inputs, d_outputs, d_residuals, mode) :

[Optional] If the user wants to implement partial derivatives in a matrix-free way, this method performs the matrix-vector product. If mode is ‘fwd’, this method computes \(d\_residuals=J⋅[d\_inputs \ d\_outputs]^T\) . If mode is ‘rev’, this method computes \([d\_inputs \ d\_outputs]^T=J^T⋅d\_residuals\).

def apply_linear(self, inputs, outputs,
                 d_inputs, d_outputs, d_residuals, mode):
    a = inputs['a']
    b = inputs['b']
    c = inputs['c']
    x = outputs['x']
    if mode == 'fwd':
        if 'x' in d_residuals:
            if 'x' in d_outputs:
                d_residuals['x'] += (2 * a * x + b) * d_outputs['x']
            if 'a' in d_inputs:
                d_residuals['x'] += x ** 2 * d_inputs['a']
            if 'b' in d_inputs:
                d_residuals['x'] += x * d_inputs['b']
            if 'c' in d_inputs:
                d_residuals['x'] += d_inputs['c']
    elif mode == 'rev':
        if 'x' in d_residuals:
            if 'x' in d_outputs:
                d_outputs['x'] += (2 * a * x + b) * d_residuals['x']
            if 'a' in d_inputs:
                d_inputs['a'] += x ** 2 * d_residuals['x']
            if 'b' in d_inputs:
                d_inputs['b'] += x * d_residuals['x']
            if 'c' in d_inputs:
                d_inputs['c'] += d_residuals['x']

solve_linear(d_outputs, d_residuals, mode) :

[Optional] Solves a linear system where the matrix is \(d\_residuals/d\_outputs\) or its transpose. If mode is ‘fwd’, the right-hand side vector is \(d\_residuals\) and the solution vector is \(d\_outputs\). If mode is ‘rev’, the right-hand side vector is \(d\_outputs\) and the solution vector is \(d\_residuals\).

def solve_linear(self, d_outputs, d_residuals, mode):
    if mode == 'fwd':
        d_outputs['x'] = self.inv_jac * d_residuals['x']
    elif mode == 'rev':
        d_residuals['x'] = self.inv_jac * d_outputs['x']

guess_nonlinear(self, inputs, outputs, residuals) :

[Optional] This method allows the user to calculate and specify an initial guess for implicit states. It is called at the start of the solve loop from certain nonlinear solvers (i.e. NewtonSolver and BroydenSolver), so it is useful for when you would like to “reset” the initial conditions on an inner-nested solve whenever an outer loop solver or driver changes other values. However if you include this function without checking the residuals, you will be starting the solvers from scratch during each run_model during an optimization. Add a check of your residuals and you can skip it when the system is converged.

Since it is a hook for custom code, you could also use it to monitor variables in the input, output, or residual vectors and change the initial guess when some condition is met.

Here is a simple example where we use NewtonSolver to find one of the roots of a second-order quadratic equation. Which root you get depends on the initial guess.

Note

When a solver is restarted from the last successful solve due to having the restart_from_sucessful option set to True, the guess_nonlinear function will not be called.

import openmdao.api as om
import numpy as np

class ImpWithInitial(om.ImplicitComponent):
    """
    An implicit component to solve the quadratic equation: x^2 - 4x + 3
    (solutions at x=1 and x=3)
    """
    def setup(self):
        self.add_input('a', val=1.)
        self.add_input('b', val=-4.)
        self.add_input('c', val=3.)

        self.add_output('x', val=0.)

        self.declare_partials(of='*', wrt='*')

    def apply_nonlinear(self, inputs, outputs, residuals):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        x = outputs['x']
        residuals['x'] = a * x ** 2 + b * x + c

    def linearize(self, inputs, outputs, partials):
        a = inputs['a']
        b = inputs['b']
        c = inputs['c']
        x = outputs['x']

        partials['x', 'a'] = x ** 2
        partials['x', 'b'] = x
        partials['x', 'c'] = 1.0
        partials['x', 'x'] = 2 * a * x + b

    def guess_nonlinear(self, inputs, outputs, resids):
        # Check residuals
        if np.abs(resids['x']) > 1.0E-2:
            # Default initial state of zero for x takes us to x=1 solution.
            # Here we set it to a value that will take us to the x=3 solution.
            outputs['x'] = 5.0

prob = om.Problem()
model = prob.model

model.add_subsystem('comp', ImpWithInitial())

model.nonlinear_solver = om.NewtonSolver(solve_subsystems=False)
model.linear_solver = om.ScipyKrylov()

prob.setup()
prob.run_model()

print(prob.get_val('comp.x'))

NL: Newton Converged in 6 iterations
[3.]

ImplicitComponent

Contents

ImplicitComponent#

ImplicitComponent Methods#