Building Models with Solvers and Implicit Components

This tutorial will show you how to define implicit components and build models with them. We’ll use a nonlinear circuit analysis example problem.

Circuit analysis

Consider a simple electrical circuit made up from two resistors, a diode, and a constant current source. Our goal is to solve for the steady-state voltages at node 1 and node 2.

diagram of a simple circuit with two resistors and one diode

In order to find the voltages, we’ll employ Kirchoff’s current law, and solve for the voltages needed at each node to drive the net current to 0.

This means that the voltages at each node are state variables for the analysis. In other words, V1 and V2 are defined implicitly by the following residual equation:

\[\mathcal{R_{node_j}} = \sum_k I_{k}^{in} - \sum_k I_{k}^{out} = 0 .\]

To build this model we’re going to define three different components:

  1. Resistor (Explicit)
  2. Diode (Explicit)
  3. Node (Implicit)

ExplicitComponents - Resistor and Diode

The Resistor and Diode components will each compute their current, given the voltages on either side. These calculations are analytic functions, so we’ll inherit from ExplicitComponent. These components will each declare some options to allow you to pass in the relevant physical constants, and to allow you to give some reasonable default values.

class Resistor(ExplicitComponent):
    """Computes current across a resistor using Ohm's law."""

    def initialize(self):
        self.options.declare('R', default=1., desc='Resistance in Ohms')

    def setup(self):
        self.add_input('V_in', units='V')
        self.add_input('V_out', units='V')
        self.add_output('I', units='A')

        self.declare_partials('I', 'V_in', method='fd')
        self.declare_partials('I', 'V_out', method='fd')

    def compute(self, inputs, outputs):
        deltaV = inputs['V_in'] - inputs['V_out']
        outputs['I'] = deltaV / self.options['R']
class Diode(ExplicitComponent):
    """Computes current across a diode using the Shockley diode equation."""

    def initialize(self):
        self.options.declare('Is', default=1e-15, desc='Saturation current in Amps')
        self.options.declare('Vt', default=.025875, desc='Thermal voltage in Volts')

    def setup(self):
        self.add_input('V_in', units='V')
        self.add_input('V_out', units='V')
        self.add_output('I', units='A')

        self.declare_partials('I', 'V_in', method='fd')
        self.declare_partials('I', 'V_out', method='fd')

    def compute(self, inputs, outputs):
        deltaV = inputs['V_in'] - inputs['V_out']
        Is = self.options['Is']
        Vt = self.options['Vt']
        outputs['I'] = Is * np.exp(deltaV / Vt - 1)

Note

Since we’ve provided default values for the options, they won’t be required arguments when instantiating Resistor or Diode. Check out the Features section for more details on how to use component options.

ImplicitComponent - Node

The Node component inherits from ImplicitComponent, which has a different interface than ExplicitComponent. Rather than compute the values of its outputs, it computes residuals via the apply_nonlinear method. When those residuals have been driven to zero, the values of the outputs will be implicitly known. apply_nonlinear computes the residuals using values from inputs and outputs. Notice that we still define V as an output of the Node component, albeit one that is implicitly defined.

class Node(ImplicitComponent):
    """Computes voltage residual across a node based on incoming and outgoing current."""

    def initialize(self):
        self.options.declare('n_in', default=1, types=int, desc='number of connections with + assumed in')
        self.options.declare('n_out', default=1, types=int, desc='number of current connections + assumed out')

    def setup(self):
        self.add_output('V', val=5., units='V')

        for i in range(self.options['n_in']):
            i_name = 'I_in:{}'.format(i)
            self.add_input(i_name, units='A')

        for i in range(self.options['n_out']):
            i_name = 'I_out:{}'.format(i)
            self.add_input(i_name, units='A')

        #note: we don't declare any partials wrt `V` here,
        #      because the residual doesn't directly depend on it
        self.declare_partials('V', 'I*', method='fd')

    def apply_nonlinear(self, inputs, outputs, residuals):
        residuals['V'] = 0.
        for i_conn in range(self.options['n_in']):
            residuals['V'] += inputs['I_in:{}'.format(i_conn)]
        for i_conn in range(self.options['n_out']):
            residuals['V'] -= inputs['I_out:{}'.format(i_conn)]

All implicit components must define the apply_nonlinear method, but it is not a requirement that every ImplicitComponent define the solve_nonlinear method. In fact, for the Node component, it is not even possible to define a solve_nonlinear because V does not show up directly in the residual function. So the implicit function represented by instances of the Node component must be converged at a higher level in the model hierarchy.

There are cases where it is possible, and even advantageous, to define the solve_nonlinear method. For example, when a component is performing an engineering analysis with its own specialized nonlinear solver routines (e.g. CFD or FEM), then it makes sense to expose those to OpenMDAO via solve_nonlinear so OpenMDAO can make use of them. Just remember that apply_nonlinear must be defined, regardless of whether you also define solve_nonlinear.

Note

In this case, the residual equation is not a direct function of the state variable V. Often, however, the residual might be a direct function of one or more output variables. If that is the case, you can access the values via outputs['V']. See the ImplicitComponent documentation for an example of this.

Building the Circuit Group and Solving It with NewtonSolver

We can combine the Resistor, Diode, and Node into the circuit pictured above using a Group. Adding components and connecting their variables is the same as what you’ve seen before in the Sellar - Two Discipline tutorial. What is new here is the additional use of the nonlinear NewtonSolver and linear DirectSolver to converge the system.

In previous tutorials, we used a gradient-free NonlinearBlockGaussSeidel solver, but that won’t work here. Just above, we discussed that the Node class does not, and in fact can not, define its own solve_nonlinear method. Hence, there would be no calculations for the GaussSeidel solver to iterate on. Instead we use the Newton solver at the Circuit level, which uses Jacobian information to compute group level updates for all the variables simultaneously. The Newton solver’s use of that Jacobian information is why we need to declare a linear solver in this case.

Note

OpenMDAO provides a library of linear solvers that are useful in different advanced scenarios. For many problems, especially problems built from components with mostly scalar variables, the DirectSolver will be both the most efficient and the easiest to use. We recommend you stick with DirectSolver unless you have a good reason to switch.

from openmdao.api import Group, NewtonSolver, DirectSolver, Problem, IndepVarComp

from openmdao.test_suite.test_examples.test_circuit_analysis import Resistor, Diode, Node

class Circuit(Group):

    def setup(self):
        self.add_subsystem('n1', Node(n_in=1, n_out=2), promotes_inputs=[('I_in:0', 'I_in')])
        self.add_subsystem('n2', Node())  # leaving defaults

        self.add_subsystem('R1', Resistor(R=100.), promotes_inputs=[('V_out', 'Vg')])
        self.add_subsystem('R2', Resistor(R=10000.))
        self.add_subsystem('D1', Diode(), promotes_inputs=[('V_out', 'Vg')])

        self.connect('n1.V', ['R1.V_in', 'R2.V_in'])
        self.connect('R1.I', 'n1.I_out:0')
        self.connect('R2.I', 'n1.I_out:1')

        self.connect('n2.V', ['R2.V_out', 'D1.V_in'])
        self.connect('R2.I', 'n2.I_in:0')
        self.connect('D1.I', 'n2.I_out:0')

        self.nonlinear_solver = NewtonSolver()
        self.nonlinear_solver.options['iprint'] = 2
        self.nonlinear_solver.options['maxiter'] = 20
        self.linear_solver = DirectSolver()

p = Problem()
model = p.model

model.add_subsystem('ground', IndepVarComp('V', 0., units='V'))
model.add_subsystem('source', IndepVarComp('I', 0.1, units='A'))
model.add_subsystem('circuit', Circuit())

model.connect('source.I', 'circuit.I_in')
model.connect('ground.V', 'circuit.Vg')

p.setup()

# set some initial guesses
p['circuit.n1.V'] = 10.
p['circuit.n2.V'] = 1.

p.run_model()
=======
circuit
=======
NL: Newton 0 ; 21.5152153 1
NL: Newton 1 ; 8.236021 0.382799841
NL: Newton 2 ; 3.02968362 0.140815863
NL: Newton 3 ; 1.11434147 0.0517931826
NL: Newton 4 ; 0.40971227 0.0190429082
NL: Newton 5 ; 0.150488221 0.00699450221
NL: Newton 6 ; 0.0551231446 0.00256205406
NL: Newton 7 ; 0.0200406862 0.000931465752
NL: Newton 8 ; 0.0071383729 0.000331782546
NL: Newton 9 ; 0.00240366903 0.000111719497
NL: Newton 10 ; 0.00069274325 3.21978303e-05
NL: Newton 11 ; 0.000129518518 6.01985693e-06
NL: Newton 12 ; 7.66202506e-06 3.56121236e-07
NL: Newton 13 ; 3.16310628e-08 1.4701718e-09
NL: Newton 14 ; 1.15199744e-12 5.35433844e-14
NL: Newton Converged
print(p['circuit.n1.V'])
[ 9.90830282]
print(p['circuit.n2.V'])
[ 0.73858486]
print(p['circuit.R1.I'])
[ 0.09908303]
print(p['circuit.R2.I'])
[ 0.00091697]
print(p['circuit.D1.I'])
[ 0.00091697]
# sanity check: should sum to .1 Amps
print(p['circuit.R1.I'] + p['circuit.D1.I'])
[ 0.1]

Modifying Solver Settings in Your Run Script

In the above run script, we set some initial guess values: prob['n1.V']=10 and prob['n2.V']=1. If you try to play around with those initial guesses a bit, you will see that convergence is really sensitive to the initial guess you used for n2.V. Below we provide a second run script that uses the same Circuit group we defined previously, but which additionally modifies some solver settings and initial guesses. If we set the initial guess for prob['n2.V']=1e-3, then the model starts out with a massive residual. It also converges much more slowly, so although we gave it more than twice the number of iterations, it doesn’t even get close to a converged answer.

from openmdao.api import Problem, IndepVarComp

from openmdao.test_suite.test_examples.test_circuit_analysis import Circuit

p = Problem()
model = p.model

model.add_subsystem('ground', IndepVarComp('V', 0., units='V'))
model.add_subsystem('source', IndepVarComp('I', 0.1, units='A'))
model.add_subsystem('circuit', Circuit())

model.connect('source.I', 'circuit.I_in')
model.connect('ground.V', 'circuit.Vg')

p.setup()

# you can change the NewtonSolver settings in circuit after setup is called
newton = p.model.circuit.nonlinear_solver
newton.options['maxiter'] = 50

# set some initial guesses
p['circuit.n1.V'] = 10.
p['circuit.n2.V'] = 1e-3

p.run_model()
=======
circuit
=======
NL: Newton 0 ; 2.53337743 1
NL: Newton 1 ; 2.5649167e+152 1.01244949e+152
NL: Newton 2 ; 9.43598353e+151 3.72466551e+151
NL: Newton 3 ; 3.47137142e+151 1.37025434e+151
NL: Newton 4 ; 1.27707085e+151 5.04098142e+150
NL: Newton 5 ; 4.6981719e+150 1.85450926e+150
NL: Newton 6 ; 1.72839425e+150 6.82249013e+149
NL: Newton 7 ; 6.35852998e+149 2.50990236e+149
NL: Newton 8 ; 2.33921766e+149 9.23359318e+148
NL: Newton 9 ; 8.60566712e+148 3.39691474e+148
NL: Newton 10 ; 3.16590919e+148 1.24967924e+148
NL: Newton 11 ; 1.16469541e+148 4.59740185e+147
NL: Newton 12 ; 4.28475775e+147 1.69132231e+147
NL: Newton 13 ; 1.57630474e+147 6.22214727e+146
NL: Newton 14 ; 5.79901314e+146 2.28904429e+146
NL: Newton 15 ; 2.13337894e+146 8.42108607e+145
NL: Newton 16 ; 7.84841417e+145 3.0980043e+145
NL: Newton 17 ; 2.88732601e+145 1.13971411e+145
NL: Newton 18 ; 1.0622084e+145 4.19285493e+144
NL: Newton 19 ; 3.90772185e+144 1.54249494e+144
NL: Newton 20 ; 1.43759831e+144 5.67463139e+143
NL: Newton 21 ; 5.28873081e+143 2.08762056e+143
NL: Newton 22 ; 1.94565293e+143 7.68007526e+142
NL: Newton 23 ; 7.15779545e+142 2.82539639e+142
NL: Newton 24 ; 2.63325667e+142 1.03942533e+142
NL: Newton 25 ; 9.68739711e+141 3.82390598e+141
NL: Newton 26 ; 3.56386309e+141 1.40676358e+141
NL: Newton 27 ; 1.3110973e+141 5.17529398e+140
NL: Newton 28 ; 4.82335061e+140 1.90392105e+140
NL: Newton 29 ; 1.77444581e+140 7.00426945e+139
NL: Newton 30 ; 6.52794748e+139 2.57677652e+139
NL: Newton 31 ; 2.40154407e+139 9.47961423e+138
NL: Newton 32 ; 8.83495762e+138 3.48742257e+138
NL: Newton 33 ; 3.25026208e+138 1.28297586e+138
NL: Newton 34 ; 1.1957277e+138 4.71989561e+137
NL: Newton 35 ; 4.39892138e+137 1.73638611e+137
NL: Newton 36 ; 1.61830401e+137 6.38793095e+136
NL: Newton 37 ; 5.95352277e+136 2.35003387e+136
NL: Newton 38 ; 2.19022095e+136 8.64545854e+135
NL: Newton 39 ; 8.05752829e+135 3.18054791e+135
NL: Newton 40 ; 2.96425628e+135 1.1700808e+135
NL: Newton 41 ; 1.09051002e+135 4.30456987e+134
NL: Newton 42 ; 4.01183967e+134 1.58359336e+134
NL: Newton 43 ; 1.47590185e+134 5.82582696e+133
NL: Newton 44 ; 5.42964441e+133 2.14324338e+133
NL: Newton 45 ; 1.99749315e+133 7.88470414e+132
NL: Newton 46 ; 7.34850863e+132 2.9006766e+132
NL: Newton 47 ; 2.70341748e+132 1.06711991e+132
NL: Newton 48 ; 9.94550931e+131 3.92579061e+131
NL: Newton 49 ; 3.65881911e+131 1.44424556e+131
NL: Newton 50 ; 1.34603034e+131 5.31318517e+130
NL: Newton Failed to Converge in 50 iterations
print(p['circuit.n1.V'])
[ 9.98744708]
print(p['circuit.n2.V'])
[ 8.73215484]
# sanity check: should sum to .1 Amps
print(p['circuit.R1.I'] + p['circuit.D1.I'])
[ 0.09987447]

Note

You actually can get this model to converge. But you have to set the options for maxiter=400 and rtol=1e-100. (The rtol value needs to be very low to prevent premature termination.)

Tweaking Newton Solver Settings to Get More Robust Convergence

The NewtonSolver has a lot of features that allow you to modify its behavior to handle more challenging problems. We’re going to look at two of the most important ones here:

  1. Line searches
  2. The solve_subsystems option

If we use both of these in combination, we can dramatically improve the robustness of the solver for this problem. The linesearch attribute makes sure that the solver doesn’t take too big of a step. The solve_subsystems option allows the Resistor and Diode components (the two ExplicitComponents) to help the convergence by updating their own output values given their inputs. When you use NewtonSolver on models with a lot of ExplicitComponents, you may find that turning on solve_subsystems helps convergence, but you need to be careful about the execution order when you try this.

Note

For this case, we used the ArmijoGoldsteinLS, which basically limits step sizes so that the residual always goes down. For many problems you might want to use BoundsEnforceLS instead, which only activates the line search to enforce upper and lower bounds on the outputs in the model.

from openmdao.api import ArmijoGoldsteinLS, Problem, IndepVarComp

from openmdao.test_suite.test_examples.test_circuit_analysis import Circuit

p = Problem()
model = p.model

model.add_subsystem('ground', IndepVarComp('V', 0., units='V'))
model.add_subsystem('source', IndepVarComp('I', 0.1, units='A'))
model.add_subsystem('circuit', Circuit())

model.connect('source.I', 'circuit.I_in')
model.connect('ground.V', 'circuit.Vg')

p.setup()

# you can change the NewtonSolver settings in circuit after setup is called
newton = p.model.circuit.nonlinear_solver
newton.options['iprint'] = 2
newton.options['maxiter'] = 10
newton.options['solve_subsystems'] = True
newton.linesearch = ArmijoGoldsteinLS()
newton.linesearch.options['maxiter'] = 10
newton.linesearch.options['iprint'] = 2

# set some initial guesses
p['circuit.n1.V'] = 10.
p['circuit.n2.V'] = 1e-3

p.run_model()
=======
circuit
=======
NL: Newton 0 ; 0.00141407214 1
|  LS: AG 0 ; 2.56438641e+152 1.81347636e+155
|  LS: AG 1 ; 3.13138636e+68 2.21444598e+71
|  LS: AG 3 ; 4.42844906e+68 3.13169953e+71
|  LS: AG 5 ; 4.89359754e+26 3.46064207e+29
|  LS: AG 7 ; 514418.461 363785161
|  LS: AG 9 ; 0.00136161968 0.962906799
NL: Newton 1 ; 0.00136988234 0.968749966
|  LS: AG 0 ; 2.53928024e+152 1.85364842e+155
|  LS: AG 1 ; 1.3058377e+71 9.53248072e+73
|  LS: AG 3 ; 1.84673338e+71 1.34809635e+74
|  LS: AG 5 ; 4.18787282e+30 3.05710403e+33
|  LS: AG 7 ; 1.99428915e+10 1.45581054e+13
|  LS: AG 9 ; 1.37554356 1004.13264
NL: Newton 2 ; 0.00132137032 0.934443361
|  LS: AG 0 ; 3.01803415e+32 2.28401841e+35
|  LS: AG 1 ; 4.9393708e+13 3.738067e+16
|  LS: AG 3 ; 6.98532517e+13 5.28642504e+16
|  LS: AG 5 ; 28258.9479 21386092.5
|  LS: AG 7 ; 0.566402375 428.647719
|  LS: AG 9 ; 0.00207446053 1.56993122
NL: Newton 3 ; 0.00121787752 0.861255579
|  LS: AG 0 ; 0.0794923223 65.2711959
|  LS: AG 1 ; 0.0026724604 2.19435892
|  LS: AG 3 ; 0.0031239501 2.56507739
|  LS: AG 5 ; 0.000928552052 0.762434675
NL: Newton 4 ; 0.00103096606 0.729075998
|  LS: AG 0 ; 0.00196602623 1.90697474
|  LS: AG 1 ; 0.000592514378 0.574717638
NL: Newton 5 ; 0.000399320609 0.282390549
|  LS: AG 0 ; 5.7621172e-07 0.00144298017
NL: Newton 6 ; 5.76211723e-07 0.000407483965
|  LS: AG 0 ; 1.90981024e-10 0.000331442448
NL: Newton 7 ; 1.90981023e-10 1.35057482e-07
|  LS: AG 0 ; 3.69994865e-15 1.93733838e-05
NL: Newton 8 ; 3.69994865e-15 2.61652043e-12
NL: Newton Converged
print(p['circuit.n1.V'])
[ 9.90830282]
print(p['circuit.n2.V'])
[ 0.73858486]
print(p['circuit.R1.I'])
[ 0.09908303]
print(p['circuit.R2.I'])
[ 0.00091697]
print(p['circuit.D1.I'])
[ 0.00091697]
# sanity check: should sum to .1 Amps
print(p['circuit.R1.I'] + p['circuit.D1.I'])
[ 0.1]

Note

This tutorial used finite difference to approximate the partial derivatives for all the components. Check out this example if you want to see the same problem solved with analytic derivatives.