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_analysiserrorFalse [True, False] ['bool'] When the option is true, a solver will reraise any AnalysisError that arises during subsolve; when false, it will continue solving.
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.
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)[source]

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.

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
=======
BROYDEN 0 ; 60.8585191 1
BROYDEN 1 ; 22.3887968 0.367882708
BROYDEN 2 ; 12.510447 0.205566076
BROYDEN 3 ; 5.98619238 0.098362439
BROYDEN 4 ; 3.04366868 0.050012204
BROYDEN 5 ; 1.51175582 0.0248404963
BROYDEN 6 ; 0.75759176 0.0124484094
BROYDEN 7 ; 0.378237612 0.00621503147
BROYDEN 8 ; 0.188999522 0.00310555572
BROYDEN 9 ; 0.0942955644 0.00154942259
BROYDEN 10 ; 0.0469601147 0.000771627627
BROYDEN 11 ; 0.0232898966 0.000382689177
BROYDEN 12 ; 0.0114571067 0.000188258059
BROYDEN 13 ; 0.00554434897 9.11022655e-05
BROYDEN 14 ; 0.00259570587 4.26514792e-05
BROYDEN 15 ; 0.0011364197 1.86731408e-05
BROYDEN 16 ; 0.000434046425 7.13205696e-06
BROYDEN 17 ; 0.000125657613 2.06474977e-06
BROYDEN 18 ; 2.10704958e-05 3.46220975e-07
BROYDEN 19 ; 1.29557101e-06 2.12882441e-08
BROYDEN 20 ; 1.45277217e-08 2.38713033e-10
BROYDENSolver '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
=======
BROYDEN 0 ; 60.8585191 1
BROYDEN 1 ; 22.3887968 0.367882708
BROYDEN 2 ; 12.510447 0.205566076
BROYDEN 3 ; 5.98619238 0.098362439
BROYDEN 4 ; 3.04366868 0.050012204
BROYDEN 5 ; 1.51175582 0.0248404963
BROYDEN 6 ; 0.75759176 0.0124484094
BROYDEN 7 ; 0.378237612 0.00621503147
BROYDEN 8 ; 0.188999522 0.00310555572
BROYDEN 9 ; 0.0942955644 0.00154942259
BROYDEN 10 ; 0.0469601147 0.000771627627
BROYDEN 11 ; 0.0232898966 0.000382689177
BROYDEN 12 ; 0.0114571067 0.000188258059
BROYDEN 13 ; 0.00554434897 9.11022655e-05
BROYDEN 14 ; 0.00259570587 4.26514792e-05
BROYDEN 15 ; 0.0011364197 1.86731408e-05
BROYDEN 16 ; 0.000434046425 7.13205696e-06
BROYDEN 17 ; 0.000125657613 2.06474977e-06
BROYDEN 18 ; 2.10704958e-05 3.46220975e-07
BROYDEN 19 ; 1.29557101e-06 2.12882441e-08
BROYDEN 20 ; 1.45277217e-08 2.38713033e-10
BROYDENSolver '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.

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()
BROYDEN 0 ; 11.2725705 1
|  LS: BCHK 0 ; 6.82339784 0.605309838
BROYDEN 1 ; 6.67731871 0.592351025
|  LS: BCHK 0 ; 1.31638163 0.197142249
BROYDEN 2 ; 0.00261878643 0.00023231493
|  LS: BCHK 0 ; 0.000516169776 0.197102661
BROYDEN 3 ; 6.33742047e-07 5.62198343e-08
|  LS: BCHK 0 ; 1.2489701e-07 0.197078624
BROYDEN 4 ; 6.21724894e-14 5.51537817e-15
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
=======
BROYDEN 0 ; 60.8585191 1
BROYDEN 1 ; 22.3887968 0.367882708
BROYDEN 2 ; 8.23629948 0.135335194
BROYDEN 3 ; 3.02978539 0.04978408
BROYDEN 4 ; 1.11437824 0.0183109654
BROYDEN 5 ; 0.409725119 0.00673242012
BROYDEN 6 ; 0.150492273 0.0024728218
BROYDEN 7 ; 0.0551239641 0.000905772353
BROYDEN 8 ; 0.0200403261 0.000329293686
BROYDEN 9 ; 0.00713760492 0.000117281936
BROYDEN 10 ; 0.00240281896 3.94820478e-05
BROYDEN 11 ; 0.000692025597 1.13710555e-05
BROYDEN 12 ; 0.000282811849 4.64703796e-06
BROYDEN 13 ; 6.22844815e-05 1.02343078e-06
BROYDEN 14 ; 7.64241042e-06 1.25576674e-07
BROYDEN 15 ; 2.45096367e-07 4.02731402e-09
BROYDEN 16 ; 1.00716012e-09 1.65492052e-11
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 and stall_tol

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.

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()
BROYDEN 0 ; 4 1
|  LS: BCHK 0 ; 4.76340215 1.50632002
BROYDEN 1 ; 0.7 0.175
|  LS: BCHK 0 ; 0.63953754 0.913625056
BROYDEN 2 ; 0.7 0.175
|  LS: BCHK 0 ; 0.63953754 0.913625056
BROYDEN 3 ; 0.7 0.175
|  LS: BCHK 0 ; 0.63953754 0.913625056
BROYDEN 4 ; 0.7 0.175
BROYDENSolver 'BROYDEN' on system '' stalled after 4 iterations.