Components with Implicit States¶
In this tutorial, we show how to build a component that contains an implicit function. So far, we have learned how to define an OpenMDAO component that represents an explicit function of unknowns with respect to its params.
y = F(x)
OpenMDAO also allows us to define a component that contains an implicit function of the same variables:
R(x, y) = 0
Here, the variable ‘x’ is a known param that is passed in to the component, but the variable ‘y’, called the state, is an unknown that needs to be solved 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.
Some equations can easily be converted from implicit to explicit, but there are cases that are difficult or impossible to represent in an explicit form, so for that reason we support implicit equations.
In the following tutorial, we will build a component to solve the following equations:
f(x,z) = xz + z - 4 = 0
y = x + 2z
The first equation is an implicit function of the param ‘z’ and the state ‘x’. This example is fairly easy to convert to explicit, and you can do this if you would like a direct comparison between implicit and explicit components. The second equation is a normal explicit one. We will include this in the same component to show how an OpenMDAO component can contain any number of implicit and explicit relationships.
There are 3 ways to solve for the state in an implicit component.
- The component can solve it (and OpenMDAO does nothing).
- OpenMDAO can solve it (and the component does nothing).
- Both OpenMDAO and the component can solve it.
Most of the time, you will be using ‘2’ or ‘3’, mainly because if you implement ‘1’, there is no reason that OpenMDAO needs to know about the state, so it might as well be an output.
We will show you how to set up and run each of these.
Implicit Component that Solves Itself¶
So let’s implement a component with param ‘x’, output ‘y’, and state ‘z’.
from __future__ import print_function
import numpy as np
from openmdao.api import Component, Group, Problem, ScipyGMRES
class SimpleImplicitComp(Component):
""" A Simple Implicit Component with an additional output equation.
f(x,z) = xz + z - 4
y = x + 2z
Sol: when x = 0.5, z = 2.666
Coupled derivs:
y = x + 8/(x+1)
dy_dx = 1 - 8/(x+1)**2 = -2.5555555555555554
z = 4/(x+1)
dz_dx = -4/(x+1)**2 = -1.7777777777777777
"""
def __init__(self):
super(SimpleImplicitComp, self).__init__()
# Params
self.add_param('x', 0.5)
# Unknowns
self.add_output('y', 0.0)
# States
self.add_state('z', 0.0)
self.maxiter = 25
self.atol = 1.0e-12
def solve_nonlinear(self, params, unknowns, resids):
""" Simple iterative solve. (Babylonian method)."""
x = params['x']
z = unknowns['z']
znew = z
itercount = 0
eps = 1.0e99
while itercount < self.maxiter and abs(eps) > self.atol:
z = znew
znew = 4.0 - x*z
eps = x*znew + znew - 4.0
itercount += 1
# Our State
unknowns['z'] = znew
# Our Output
unknowns['y'] = x + 2.0*znew
def apply_nonlinear(self, params, unknowns, resids):
""" Don't solve; just calculate the residual."""
x = params['x']
z = unknowns['z']
resids['z'] = x*z + z - 4.0
# Output equations need to evaluate a residual just like an explicit comp.
resids['y'] = x + 2.0*z - unknowns['y']
def linearize(self, params, unknowns, resids):
"""Analytical derivatives."""
J = {}
# Output equation
J[('y', 'x')] = np.array([1.0])
J[('y', 'z')] = np.array([2.0])
# State equation
J[('z', 'z')] = np.array([params['x'] + 1.0])
J[('z', 'x')] = np.array([unknowns['z']])
return J
Since we are solving the implicit equation in our component, we include code in solve_nonlinear that iterates on the implicit equation using the Babylonian method, which is essentially fixed point iteration. The ‘solve_nonlinear’ method must also set values for the unknowns, in this case ‘y’.
When you have explicit equations, you must also specify a new method called ‘apply_nonlinear’. This method is called when OpenMDAO wants to evaluate the residuals, so in our case we want to return the current value of:
R = xz + z - 4
This is an in-place evaluation using current values of all states, params, and unknowns as they appear in the vectors. You should never set anything in this method except residuals.
The value of a residual is placed in the resids vector in the variable named for the appropriate state. In our case, the residual for state ‘x’ is placed in resids[‘x’].
If your component has unknowns, then there is one less obvious thing you need to do. For each unknown, you need to define a residual. This is done by rearranging the equation so that it is in implicit form.
y = x + 2.0*z
R = (x + 2.0*z) - y
By convention, OpenMDAO expects the current value of the output to be subtracted as shown. Note that the residuals are important so that this component can be correctly converged by solvers in the containing group; they don’t impact the self-solve case.
Finally, we show how to declare derivatives for the implicit comp. The derivatives for the output equation are as expected. For the implicit equation, we evaluate the derivatives of the state equation with respect to all inputs and states. All derivatives are assigned to the state output in the Jacobian, so the derivative of the residual with respect to the state resides in the (‘z’, ‘z’) key.
Now, let’s put the implicit component into a simple model and run it.
top = Problem()
root = top.root = Group()
root.add('comp', SimpleImplicitComp())
root.ln_solver = ScipyGMRES()
top.setup()
top.run()
print('Solution: x = %f, z = %f, y = %f' % (top['comp.x'], top['comp.z'], top['comp.y']))
Note that we need to specify ScipyGMRES as our linear solver as we need one that can handle implicit states. We aren’t actually calculating any derivatives here, but if we wanted to, for example, place this in a larger model and optimize it, GMRES would be needed here so we add it.
Solution: x = 0.500000, z = 2.666667, y = 5.833333
This matches the expected answer.
Implicit Component that is Solved by OpenMDAO¶
Coming soon.
Implicit Component that is Solved by both itself and OpenMDAO¶
Coming soon.
from __future__ import print_function
import numpy as np
from openmdao.api import Component, Group, Problem, ScipyGMRES
class SimpleImplicitComp(Component):
""" A Simple Implicit Component with an additional output equation.
f(x,z) = xz + z - 4
y = x + 2z
Sol: when x = 0.5, z = 2.666
Coupled derivs:
y = x + 8/(x+1)
dy_dx = 1 - 8/(x+1)**2 = -2.5555555555555554
z = 4/(x+1)
dz_dx = -4/(x+1)**2 = -1.7777777777777777
"""
def __init__(self):
super(SimpleImplicitComp, self).__init__()
# Params
self.add_param('x', 0.5)
# Unknowns
self.add_output('y', 0.0)
# States
self.add_state('z', 0.0)
self.maxiter = 25
self.atol = 1.0e-3
def solve_nonlinear(self, params, unknowns, resids):
""" Simple iterative solve. (Babylonian method)."""
x = params['x']
z = unknowns['z']
znew = z
itercount = 0
eps = 1.0e99
while itercount < self.maxiter and abs(eps) > self.atol:
z = znew
znew = 4.0 - x*z
eps = x*znew + znew - 4.0
itercount += 1
# Our State
unknowns['z'] = znew
# Our Output
unknowns['y'] = x + 2.0*znew
def apply_nonlinear(self, params, unknowns, resids):
""" Don't solve; just calculate the residual."""
x = params['x']
z = unknowns['z']
resids['z'] = x*z + z - 4.0
# Output equations need to evaluate a residual just like an explicit comp.
resids['y'] = x + 2.0*z - unknowns['y']
def linearize(self, params, unknowns, resids):
"""Analytical derivatives."""
J = {}
# Output equation
J[('y', 'x')] = np.array([1.0])
J[('y', 'z')] = np.array([2.0])
# State equation
J[('z', 'z')] = np.array([params['x'] + 1.0])
J[('z', 'x')] = np.array([unknowns['z']])
return J