# Implicit Components¶

In the basic tutorial, we learned how to define an OpenMDAO component that represents an explicit function of outputs with respect to its inputs.

$y = F(x)$

OpenMDAO also allows you to define a component that represents an implicit function of the same variables:

$R(x, y) = 0$

Here, the variable x is a known input that is passed to the implicit component, but the variable y, called the state, is an unknown that needs to be solved in order to satisfy the equation. The left hand side of the implicit equation is called the residual. Since an implicit function may not have a closed-form solution, the state is typically determined by numerically solving the residual equation, or in other words, iterating on the state until the residual is driven to zero.

An implicit equation can be solved using the base Component class in OpenMDAO as long as the component solves the residual equations in its execute method. However, there are some advantages to letting OpenMDAO be aware of states and residuals. First, this allows you to define analytic derivatives so that the implicit component can form part of the coupled solution of the full model gradient without having to resort to finite difference. In addition, solving the residual equations can be performed by OpenMDAO, both locally as part of the execution of the implicit component, as well as globally if you want to converge your implicit component simultaneously with other solver loops in the model.

For the sake of convenience, we’ve broadened the definition of an implicit component to allow the inclusion of additional inputs and outputs that aren’t part of the implicit equations.

$R(x, y) = 0$$z = F(x, y, u)$

In this set of equations, u is an input that does not affect the residuals, and ‘z’ is an output that is not a state, but satisfies an additional explicit equation. This is analogous to the output equation of a state space model in control theory.

Let’s set up a simple component that can solve this set of equations, three of which are implicit while one is explicit:

$r_0(x, y, z) = c*(3x + 2y - z) - 1 = 0$$r_1(x, y, z) = 2x - 2y + 4z + 2 = 0$$r_2(x, y, z) = -x + y/2. - z =0$$y_{out} = c + x + y + z$

In these equations, the states are x, y, and z, and the residuals are r_0, r_1, and r_2. The variable c is a normal input, and y_out is an explicit output. Note that the number of states must equal the number of residuals in order for the system to have a unique and valid solution. You can start defining your implicit component by inheriting from ImplicitComponent instead of Component. This class allows the definition of two additional iotypes: state and residual.

When an ImplicitComponent executes, it must solve its residual equations to find its states. You can provide this method, or you can use the one built into the ImplicitComponent base class, which calls scipy.optimize.fsolve. We will use the built-in here, and don’t need to define a new method called evaluate which is called during the iterative solution. The evaluate function needs to calculate the residuals for the current value of the state. It should also calculate any explicit outputs.

Let’s write our first implicit component.

import numpy as np

from openmdao.main.api import ImplicitComponent
from openmdao.main.datatypes.api import Float, Array

class MyComp_No_Deriv(ImplicitComponent):
''' Single implicit component with 3 states and residuals.

For c=2.0, (x,y,z) = (1.0, -2.333333, -2.1666667)
'''

# External inputs
c = Float(2.0, iotype="in", fd_step = .00001,
desc="non-state input")

# States
x = Float(0.0, iotype="state")
y = Float(0.0, iotype="state")
z = Float(0.0, iotype="state")

# Residuals
res = Array(np.zeros((3)), iotype="residual")

# Outputs
y_out = Float(iotype='out')

def evaluate(self):
"""run a single step to calculate the residual
values for the given state var values"""

c, x, y, z = self.c, self.x, self.y, self.z

self.res = self.c*(3*x + 2*y - z) - 1
self.res = 2*x - 2*y + 4*z + 2
self.res = -x + y/2. - z

self.y_out = c + x + y + z


We have taken our 3 residuals and placed them in a single variable array called res, but you could also create a separate floating point variable for each of them. Also, the initial value of our states serves as their initial conditions for their iterative solution. Now, let’s put this in an assembly:

from openmdao.main.api import Assembly, set_as_top

class Model(Assembly):

def configure(self):


and run the model. We will let the implicit component solve its own residuals.

>>> top = set_as_top(Model())
>>> top.run()
>>> # The residuals will vary depending on your system, but should be near zero.
>>> print top.comp.res
[...]
>>> print top.comp.x, top.comp.y, top.comp.z
1.0 -2.3333... -2.1666...


The implicit component completes its iteration until the state values satisfy the residual equations. We can also configure an OpenMDAO solver to solve for the states. Here, we set up a new assembly with the Brent solver as the top driver. Then we assign the states as the solver’s parameters, and constrain the residuals to be equal to zero. Also, we don’t want the implicit component’s internal solver to solve this in competition with the BroydenSolver solver, so we set eval_only to True. This means that running the implicit component just runs the eval statement we defined in the class definition.

from openmdao.main.api import Assembly, set_as_top
from openmdao.lib.drivers.api import BroydenSolver

class Model2(Assembly):

def configure(self):
self.comp.eval_only = True


Now, when we run the model, we get the same solution for the state.

>>> top = set_as_top(Model2())
>>> top.run()
>>> # The residuals will vary depending on your system, but should be near zero.
>>> print top.comp.res
[...]
>>> print top.comp.x, top.comp.y, top.comp.z
1.0 -2.3333... -2.1666...


Finally, since one of the advantageous to this implementation of implicit components is in the derivative calculation, lets specify the analytic derivatives for this simple set of equations using the apply_deriv and apply_derivT methods. To do this, we need to provide all permutations of the derivatives: namely, the derivatives of the residuals with respect to both the states and the explicit inputs, and the derivatives of the explicit output with respect to both the states and the explicit inputs. Here, we specify these as separate Jacobians in the linearize method, but this was purely to make the matrix-vector multiplication in apply_deriv and apply_derivT clean and simple.

class MyComp_Deriv(MyComp_No_Deriv):
''' This time with derivatives.
'''

def linearize(self):
#partial w.r.t c
c, x, y, z = self.c, self.x, self.y, self.z

dc = [3*x + 2*y - z, 0, 0]
dx = [3*c, 2, -1]
dy = [2*c, -2, .5]
dz = [-c, 4, -1]

self.J_res_state = np.array([dx, dy, dz]).T
self.J_res_input = np.array([dc]).T

self.J_output_input = np.array([[1.0]])
self.J_output_state = np.array([[1.0, 1.0, 1.0]])

def apply_deriv(self, arg, result):

# Residual Equation derivatives
res = self.list_residuals()
if res in result:

# wrt States
for k, state in enumerate(self.list_states()):
if state in arg:
result[res] += self.J_res_state[:, k]*arg[state]

# wrt External inputs
for k, inp in enumerate(['c']):
if inp in arg:
result[res] += self.J_res_input[:, k]*arg[inp]

# Output Equation derivatives
for j, outp in enumerate(['y_out']):
if outp in result:

# wrt States
for k, state in enumerate(self.list_states()):
if state in arg:
result[outp] += self.J_output_state[j, k]*arg[state]

# wrt External inputs
for k, inp in enumerate(['c']):
if inp in arg:
result[outp] += self.J_output_input[j, k]*arg[inp]

def apply_derivT(self, arg, result):

# wrt States
for k, state in enumerate(self.list_states()):
if state in result:

# Residual Equation derivatives
res = self.list_residuals()
if res in arg:
result[state] += self.J_res_state.T[k, :].dot(arg[res])

# Output Equation derivatives
for j, outp in enumerate(['y_out']):
if outp in arg:
result[state] += self.J_output_state.T[k, j]*arg[outp]

# wrt External inputs
for k, inp in enumerate(['c']):
if inp in result:

# Residual Equation derivatives
res = self.list_residuals()
if res in arg:
result[inp] += self.J_res_input.T[k, :].dot(arg[res])

# Output Equation derivatives
for j, outp in enumerate(['y_out']):
if outp in arg:
result[inp] += self.J_output_input.T[k, j]*arg[outp]


Specifying these derivative function removes the need for finite differencing this component in any workflow. #### Previous topic

Summary: Using the Iteration Hierarchy

#### Next topic

Working with Parallel Computational Resources