BroydenSolver

BroydenSolver is a quasi-Newton solver that implements Broyden’s second method to solve for values of the model’s states that drive their residuals to zero. It does so by maintaining an approximation to the inverse of the Jacobian of the model or a subset of the model. In some cases this can be more efficient than NewtonSolver because updating the approximated inverse Jacobian is cheaper than solving the linear system. It may take more iterations because the search direction depends on an approximation, but the iterations take fewer operations.

The BroydenSolver has two different modes of operation. It can operate on the entire model and solve for every state in the containing system and all subsystems. Alternatively, it can operate on a subset of the model and only solve for a list of states that you provide. The advantage of full-model mode is that you don’t have to worry about forgetting a state, particularly in large models where you might not be familiar with every component or variable. The disadvantage is that you are computing the inverse of a larger matrix every time you recalculate the inverse jacobian, though ideally you are not recomputing this very often. Operating on a subset of states is more efficient in both the linear solve and the Broyden update, but you do run the risk of missing a state. The BroydenSolver will print a warning if it finds any states in the model that aren’t covered by a solver.

BroydenSolver Options

Option	Default	Acceptable Values	Acceptable Types	Description
alpha	0.4	N/A	N/A	Value to scale the starting Jacobian, which is Identity. This option does nothing if you compute the initial Jacobian instead.
atol	1e-10	N/A	N/A	absolute error tolerance
compute_jacobian	True	[True, False]	['bool']	When True, compute an initial Jacobian, otherwise start with Identity scaled by alpha. Further Jacobians may also be computed depending on the other options.
converge_limit	1.0	N/A	N/A	Ratio of current residual to previous residual above which the convergence is considered a failure. The Jacobian will be regenerated once this condition has been reached a number of consecutive times as specified in max_converge_failures.
cs_reconverge	True	[True, False]	['bool']	When True, when this driver solves under a complex step, nudge the Solution vector by a small amount so that it reconverges.
debug_print	False	[True, False]	['bool']	If true, the values of input and output variables at the start of iteration are printed and written to a file after a failure to converge.
diverge_limit	2.0	N/A	N/A	Ratio of current residual to previous residual above which the Jacobian will be immediately regenerated.
err_on_non_converge	False	[True, False]	['bool']	When True, AnalysisError will be raised if we don't converge.
iprint	1	N/A	['int']	whether to print output
max_converge_failures	3	N/A	N/A	The number of convergence failures before regenerating the Jacobian.
max_jacobians	10	N/A	N/A	Maximum number of jacobians to compute.
maxiter	10	N/A	['int']	maximum number of iterations
reraise_child_analysiserror	False	[True, False]	['bool']	When the option is true, a solver will reraise any AnalysisError that arises during subsolve; when false, it will continue solving.
restart_from_successful	False	[True, False]	['bool']	If True, the states are cached after a successful solve and used to restart the solver in the case of a failed solve.
rtol	1e-10	N/A	N/A	relative error tolerance
stall_limit	0	N/A	N/A	Number of iterations after which, if the residual norms are identical within the stall_tol, then terminate as if max iterations were reached. Default is 0, which disables this feature.
stall_tol	1e-12	N/A	N/A	When stall checking is enabled, the threshold below which the residual norm is considered unchanged.
stall_tol_type	rel	['abs', 'rel']	N/A	Specifies whether the absolute or relative norm of the residual is used for stall detection.
state_vars	[]	N/A	N/A	List of the state-variable/residuals that are to be solved here.
update_broyden	True	N/A	N/A	Flag controls whether to perform Broyden update to the Jacobian. There are some applications where it may be useful to turn this off.

The BroydenSolver also contains a slot for a linear solver and a slot for a linesearch. See the linesearch section for more about these.

BroydenSolver Constructor

The call signature for the BroydenSolver constructor is:

BroydenSolver.__init__(**kwargs)

Initialize all attributes.

BroydenSolver on a Full Model

Here we show an example that uses the electrical circuit model from the advanced guide. We have replaced the NewtonSolver with a BroydenSolver, and set the maximum number of iterations to 20. We also assign a DirectSolver into the “linear_solver” slot on the BroydenSolver. This is the linear solver that will be used to assemble the Jacobian and compute its inverse. Since we don’t specify any states in the state_vars option, the BroydenSolver operates on the entire model. If you don’t specify a linear_solver here, then the BroydenSolver will use the one from the system.

Note

In this mode, only the DirectSolver can be used as the linear_solver.

Depending on the values of some of the other options such as “converge_limit”, “diverge_limit”, and “max_converge_failures”, the Jacobian might be recalculated if convergence stalls, though this doesn’t happen in the electrical circuit example.

Circuit class definition

class Circuit(om.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 = om.NewtonSolver(solve_subsystems=False)
        self.nonlinear_solver.options['iprint'] = 2
        self.nonlinear_solver.options['maxiter'] = 20
        self.linear_solver = om.DirectSolver()

import openmdao.api as om
from openmdao.test_suite.scripts.circuit_analysis import Circuit

p = om.Problem()
model = p.model

model.add_subsystem('circuit', Circuit(), promotes_inputs=[('Vg', 'V'), ('I_in', 'I')])
model.set_input_defaults('V', 0., units='V')
model.set_input_defaults('I', 0.1, units='A')

p.setup()

# Replace existing solver with BroydenSolver
model.circuit.nonlinear_solver = om.BroydenSolver()
model.circuit.nonlinear_solver.options['maxiter'] = 20
model.circuit.nonlinear_solver.linear_solver = om.DirectSolver()

# set some initial guesses
p.set_val('circuit.n1.V', 10.)
p.set_val('circuit.n2.V', 1.)

p.set_solver_print(level=2)
p.run_model()

=======
circuit
=======
NL: BROYDEN 0 ; 60.8585191 1
NL: BROYDEN 1 ; 22.3883642 0.367875599
NL: BROYDEN 2 ; 12.5102871 0.205563449
NL: BROYDEN 3 ; 5.98610005 0.0983609219
NL: BROYDEN 4 ; 3.04362478 0.0500114827
NL: BROYDEN 5 ; 1.51173342 0.0248401283
NL: BROYDEN 6 ; 0.757580648 0.0124482268
NL: BROYDEN 7 ; 0.378232039 0.0062149399
NL: BROYDEN 8 ; 0.188996739 0.00310550999
NL: BROYDEN 9 ; 0.0942941723 0.00154939972
NL: BROYDEN 10 ; 0.0469594187 0.000771616191
NL: BROYDEN 11 ; 0.0232895486 0.000382683459
NL: BROYDEN 12 ; 0.0114569328 0.000188255201
NL: BROYDEN 13 ; 0.00554426212 9.11008384e-05
NL: BROYDEN 14 ; 0.00259566267 4.26507694e-05
NL: BROYDEN 15 ; 0.00113639853 1.8672793e-05
NL: BROYDEN 16 ; 0.000434036582 7.13189523e-06
NL: BROYDEN 17 ; 0.000125653739 2.0646861e-06
NL: BROYDEN 18 ; 2.10695157e-05 3.46204869e-07
NL: BROYDEN 19 ; 1.29547486e-06 2.12866642e-08
NL: BROYDEN 20 ; 1.45259789e-08 2.38684396e-10
NL: BROYDENSolver 'NL: BROYDEN' on system 'circuit' failed to converge in 20 iterations.

print(p.get_val('circuit.n1.V'))
print(p.get_val('circuit.n2.V'))

# sanity check: should sum to .1 Amps
print(p.get_val('circuit.R1.I') + p.get_val('circuit.D1.I'))

[9.90804735]
[0.71278226]
[0.10000001]

BroydenSolver on a Subset of States

The BroydenSolver can also be used to solve for specific states. Here we consider the same circuit example, but instead we specify the two voltages n1.V’ and ‘n2.V’ as our “state_vars”. In this mode, we aren’t limited to just using the DirectSolver, and in this example we choose LinearBlockGS instead.

from openmdao.test_suite.scripts.circuit_analysis import Circuit

p = om.Problem()
model = p.model

model.add_subsystem('circuit', Circuit(), promotes_inputs=[('Vg', 'V'), ('I_in', 'I')])
model.set_input_defaults('V', 0., units='V')
model.set_input_defaults('I', 0.1, units='A')

p.setup()

# Replace existing solver with om.BroydenSolver
model.circuit.nonlinear_solver = om.BroydenSolver()
model.circuit.nonlinear_solver.options['maxiter'] = 20

# Specify states for Broyden to solve
model.circuit.nonlinear_solver.options['state_vars'] = ['n1.V', 'n2.V']

model.nonlinear_solver.linear_solver = om.LinearBlockGS()

# set some initial guesses
p.set_val('circuit.n1.V', 10.)
p.set_val('circuit.n2.V', 1.)

p.set_solver_print(level=2)
p.run_model()

=======
circuit
=======
NL: BROYDEN 0 ; 60.8585191 1
NL: BROYDEN 1 ; 22.3883642 0.367875599
NL: BROYDEN 2 ; 12.5102871 0.205563449
NL: BROYDEN 3 ; 5.98610005 0.0983609219
NL: BROYDEN 4 ; 3.04362478 0.0500114827
NL: BROYDEN 5 ; 1.51173342 0.0248401283
NL: BROYDEN 6 ; 0.757580648 0.0124482268
NL: BROYDEN 7 ; 0.378232039 0.0062149399
NL: BROYDEN 8 ; 0.188996739 0.00310550999
NL: BROYDEN 9 ; 0.0942941723 0.00154939972
NL: BROYDEN 10 ; 0.0469594187 0.000771616191
NL: BROYDEN 11 ; 0.0232895486 0.000382683459
NL: BROYDEN 12 ; 0.0114569328 0.000188255201
NL: BROYDEN 13 ; 0.00554426212 9.11008384e-05
NL: BROYDEN 14 ; 0.00259566267 4.26507694e-05
NL: BROYDEN 15 ; 0.00113639853 1.8672793e-05
NL: BROYDEN 16 ; 0.000434036582 7.13189523e-06
NL: BROYDEN 17 ; 0.000125653739 2.0646861e-06
NL: BROYDEN 18 ; 2.10695157e-05 3.46204869e-07
NL: BROYDEN 19 ; 1.29547486e-06 2.12866642e-08
NL: BROYDEN 20 ; 1.45259789e-08 2.38684396e-10
NL: BROYDENSolver 'NL: BROYDEN' on system 'circuit' failed to converge in 20 iterations.

print(p.get_val('circuit.n1.V'))
print(p.get_val('circuit.n2.V'))

[9.90804735]
[0.71278226]

# sanity check: should sum to .1 Amps
print(p.get_val('circuit.R1.I') + p.get_val('circuit.D1.I'))

[0.10000001]

BroydenSolver for Models Without Derivatives

The BroydenSolver can be used for models where you don’t have any partial derivatives defined, and don’t wish to use finite difference to calculate them. This behavior is activated by setting the “compute_jacobian” option to False. Instead of calculating an initial Jacobian, we start with an estimate that is just the identity matrix scaled by a tunable parameter in the options called “alpha”. As the BroydenSolver iterates, this estimate of the Jacobian is improved, and for some problems, a solution can be reached that satisfies the residual equations.

In this example, we solve for the coupling variable in a version of the Sellar model that severs the cycle and expresses the difference across the broken cycle as an implicit state, which the BroydenSolver will solve.

SellarStateConnection class definition

class SellarStateConnection(om.Group):
    """
    Group containing the Sellar MDA. This version uses the disciplines with derivatives.
    """

    def initialize(self):
        self.options.declare('nonlinear_solver', default=om.NewtonSolver(solve_subsystems=False),
                             desc='Nonlinear solver (class or instance) for Sellar MDA')
        self.options.declare('nl_atol', default=None,
                             desc='User-specified atol for nonlinear solver.')
        self.options.declare('nl_maxiter', default=None,
                             desc='Iteration limit for nonlinear solver.')
        self.options.declare('linear_solver', default=om.ScipyKrylov,
                             desc='Linear solver (class or instance)')
        self.options.declare('ln_atol', default=None,
                             desc='User-specified atol for linear solver.')
        self.options.declare('ln_maxiter', default=None,
                             desc='Iteration limit for linear solver.')

    def setup(self):
        sub = self.add_subsystem('sub', om.Group(),
                                 promotes=['x', 'z', 'y1',
                                           'state_eq.y2_actual', 'state_eq.y2_command',
                                           'd1.y2', 'd2.y2'])

        subgrp = sub.add_subsystem('state_eq_group', om.Group(),
                                   promotes=['state_eq.y2_actual', 'state_eq.y2_command'])
        subgrp.add_subsystem('state_eq', StateConnection())

        sub.add_subsystem('d1', SellarDis1withDerivatives(), promotes=['x', 'z', 'y1'])
        sub.add_subsystem('d2', SellarDis2withDerivatives(), promotes=['z', 'y1'])

        self.connect('state_eq.y2_command', 'd1.y2')
        self.connect('d2.y2', 'state_eq.y2_actual')

        self.add_subsystem('obj_cmp', om.ExecComp('obj = x**2 + z[1] + y1 + exp(-y2)',
                                                  z=np.array([0.0, 0.0]), x=0.0, y1=0.0, y2=0.0),
                           promotes=['x', 'z', 'y1', 'obj'])
        self.connect('d2.y2', 'obj_cmp.y2')

        self.add_subsystem('con_cmp1', om.ExecComp('con1 = 3.16 - y1'), promotes=['con1', 'y1'])
        self.add_subsystem('con_cmp2', om.ExecComp('con2 = y2 - 24.0'), promotes=['con2'])
        self.connect('d2.y2', 'con_cmp2.y2')

        self.set_input_defaults('x', 1.0)
        self.set_input_defaults('z', np.array([5.0, 2.0]))

        nl = self.options['nonlinear_solver']
        self.nonlinear_solver = nl() if inspect.isclass(nl) else nl
        if self.options['nl_atol']:
            self.nonlinear_solver.options['atol'] = self.options['nl_atol']
        if self.options['nl_maxiter']:
            self.nonlinear_solver.options['maxiter'] = self.options['nl_maxiter']

        ln = self.options['linear_solver']
        self.linear_solver = ln() if inspect.isclass(ln) else ln
        if self.options['ln_atol']:
            self.linear_solver.options['atol'] = self.options['ln_atol']
        if self.options['ln_maxiter']:
            self.linear_solver.options['maxiter'] = self.options['ln_maxiter']

    def configure(self):
        self.sub.linear_solver = om.ScipyKrylov()
        self.sub.state_eq_group.linear_solver = om.ScipyKrylov()

from openmdao.test_suite.components.sellar import SellarStateConnection

prob = om.Problem()
model = prob.model = SellarStateConnection(nonlinear_solver=om.BroydenSolver(),
                                           linear_solver=om.LinearRunOnce())

prob.setup()

model.nonlinear_solver.options['state_vars'] = ['state_eq.y2_command']
model.nonlinear_solver.options['compute_jacobian'] = False

prob.set_solver_print(level=2)
prob.run_model()

NL: BROYDEN 0 ; 11.2725705 1
|  LS: BCHK 0 ; 6.82339784 0.605309838
NL: BROYDEN 1 ; 6.67731871 0.592351025
|  LS: BCHK 0 ; 1.31638163 0.197142249
NL: BROYDEN 2 ; 0.00261878643 0.00023231493
|  LS: BCHK 0 ; 0.000516169776 0.197102661
NL: BROYDEN 3 ; 6.33742047e-07 5.62198343e-08
|  LS: BCHK 0 ; 1.2489701e-07 0.197078624
NL: BROYDEN 4 ; 6.21724894e-14 5.51537817e-15
NL: BROYDEN Converged

print(prob['y1'])
print(prob['state_eq.y2_command'])

[25.58830237]
[12.05848815]

BroydenSolver Option Examples

There are a few additional options that give you more control over when and how often the Jacobian is recomputed. The “diverge_limit” option allows you to define a limit to the ratio of current residual and the previous iteration’s residual above which the solution is considered to be diverging. If this limit is exceeded, then the Jacobian is always recomputed on the next iteration. There is also a “converge_limit” that allows you similarly define a limit above which the solution is considered to be non-converging. When this limit is exceeded, the Jacobian is not immediately recomputed until the limit has been exceeded a number of consecutive times as defined by the “max_converge_failures” option. The default value for “max_converge_failures” is 3, and the default “converge_limit” is 1.0. Exploring these options can help you solve more quickly (or in some cases solve at all) some tougher problems.

Here, we take the same circuit example from above and specify a much lower “converge_limit” and “max_converge_failures” to force recomputation of the Jacobian much more frequently. This results in a quicker convergence in terms of the number of iterations, though keep in mind that solving for the derivatives adds computational cost.

from openmdao.test_suite.scripts.circuit_analysis import Circuit

p = om.Problem()
model = p.model

model.add_subsystem('circuit', Circuit(), promotes_inputs=[('Vg', 'V'), ('I_in', 'I')])
model.set_input_defaults('V', 0., units='V')
model.set_input_defaults('I', 0.1, units='A')

p.setup()

# Replace existing solver with BroydenSolver
model.circuit.nonlinear_solver = om.BroydenSolver()
model.circuit.nonlinear_solver.options['maxiter'] = 20
model.circuit.nonlinear_solver.options['converge_limit'] = 0.1
model.circuit.nonlinear_solver.options['max_converge_failures'] = 1

# Specify states for Broyden to solve
model.circuit.nonlinear_solver.options['state_vars'] = ['n1.V', 'n2.V']

# set some initial guesses
p.set_val('circuit.n1.V', 10.)
p.set_val('circuit.n2.V', 1.)

p.set_solver_print(level=2)
p.run_model()

=======
circuit
=======
NL: BROYDEN 0 ; 60.8585191 1
NL: BROYDEN 1 ; 22.3883642 0.367875599
NL: BROYDEN 2 ; 8.23598116 0.135329963
NL: BROYDEN 3 ; 3.02960973 0.0497811937
NL: BROYDEN 4 ; 1.11429207 0.0183095495
NL: BROYDEN 5 ; 0.409685491 0.00673176898
NL: BROYDEN 6 ; 0.150474776 0.0024725343
NL: BROYDEN 7 ; 0.055116451 0.000905648902
NL: BROYDEN 8 ; 0.0200371642 0.000329241731
NL: BROYDEN 9 ; 0.00713629394 0.000117260394
NL: BROYDEN 10 ; 0.00240228312 3.94732431e-05
NL: BROYDEN 11 ; 0.000691815028 1.13675955e-05
NL: BROYDEN 12 ; 0.000282709816 4.6453614e-06
NL: BROYDEN 13 ; 6.22512906e-05 1.0228854e-06
NL: BROYDEN 14 ; 7.63617209e-06 1.25474169e-07
NL: BROYDEN 15 ; 2.44772377e-07 4.02199035e-09
NL: BROYDEN 16 ; 1.00501245e-09 1.65139156e-11
NL: BROYDEN Converged

print(p.get_val('circuit.n1.V'))
print(p.get_val('circuit.n2.V'))

[9.90804735]
[0.71278188]

# sanity check: should sum to .1 Amps
print(p.get_val('circuit.R1.I') + p.get_val('circuit.D1.I'))

[0.1]

**stall_limit, stall_tol, and stall_tol_type **

In some cases, nonlinear solvers can stall out where the norm of the residual stops changing at all. This can happen for a couple of reasons. You can hit numerical noise problems and just be wandering around in a circle, or you can get stuck on a bound and the line search just keeps running into the same spot no matter what. Either way, if you have say 100 max iterations and you stall at 15 … you waste a lot of compute time. To remedy this, you can turn on stall detection in all nonlinear solvers by setting the “stall_limit” option to a number greater than zero.

In this example, we set stall_limit to 3. While the solver iterates, it will compare the value of the residual norm to the value computed in the previous iteration. If the value matches for three iterations in a row, then iteration will terminate due to detection of a stall. If “err_on_non_converge” is set to True, then an AnalysisError will be raised just as if we had reached the iteration count limit.

We also set the stall_tol to 1e-6, which is the threshold below which a change in the relative residual norm is considered to be unchanged. The option also exists to use the absolute residual norm by setting stall_tol_type to 'abs'.

prob = om.Problem()

prob.model.add_subsystem('comp', om.ExecComp('y=3*x+1'), promotes=['*'])

balance = prob.model.add_subsystem('balance', om.BalanceComp(),
                                   promotes=['*'])
balance.add_balance('x', lower=-.1, upper=10, rhs_val=0, lhs_name='y')

nl_solver = prob.model.nonlinear_solver = om.BroydenSolver()
nl_solver.options['stall_limit'] = 3
nl_solver.options['stall_tol'] = 1e-8
nl_solver.options['maxiter'] = 100

prob.model.linear_solver = om.DirectSolver()

prob.setup()
prob.set_solver_print()

prob.run_model()

NL: BROYDEN 0 ; 4 1
|  LS: BCHK 0 ; 4.76340215 1.50632002
NL: BROYDEN 1 ; 0.7 0.175
|  LS: BCHK 0 ; 0.63953754 0.913625056
NL: BROYDEN 2 ; 0.7 0.175
|  LS: BCHK 0 ; 0.63953754 0.913625056
NL: BROYDEN 3 ; 0.7 0.175
|  LS: BCHK 0 ; 0.63953754 0.913625056
NL: BROYDEN 4 ; 0.7 0.175
NL: BROYDENSolver 'NL: BROYDEN' on system '' stalled after 4 iterations.