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¶

OptionDefaultAcceptable ValuesAcceptable TypesDescription
alpha0.4N/AN/AValue to scale the starting Jacobian, which is Identity. This option does nothing if you compute the initial Jacobian instead.
atol1e-10N/AN/Aabsolute error tolerance
compute_jacobianTrue[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_limit1.0N/AN/ARatio 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_reconvergeTrue[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_printFalse[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_limit2.0N/AN/ARatio of current residual to previous residual above which the Jacobian will be immediately regenerated.
err_on_non_convergeFalse[True, False]['bool']When True, AnalysisError will be raised if we don't converge.
iprint1N/A['int']whether to print output
max_converge_failures3N/AN/AThe number of convergence failures before regenerating the Jacobian.
max_jacobians10N/AN/AMaximum number of jacobians to compute.
maxiter10N/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.
restart_from_successfulFalse[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.
rtol1e-10N/AN/Arelative error tolerance
stall_limit0N/AN/ANumber 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_tol1e-12N/AN/AWhen stall checking is enabled, the threshold below which the residual norm is considered unchanged.
state_vars[]N/AN/AList of the state-variable/residuals that are to be solved here.
update_broydenTrueN/AN/AFlag 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
=======
NL: BROYDEN 0 ; 60.8585191 1
NL: BROYDEN 1 ; 22.3887968 0.367882708
NL: BROYDEN 2 ; 12.510447 0.205566076
NL: BROYDEN 3 ; 5.98619238 0.098362439
NL: BROYDEN 4 ; 3.04366868 0.050012204
NL: BROYDEN 5 ; 1.51175582 0.0248404963
NL: BROYDEN 6 ; 0.75759176 0.0124484094
NL: BROYDEN 7 ; 0.378237612 0.00621503147
NL: BROYDEN 8 ; 0.188999522 0.00310555572
NL: BROYDEN 9 ; 0.0942955644 0.00154942259
NL: BROYDEN 10 ; 0.0469601147 0.000771627627
NL: BROYDEN 11 ; 0.0232898966 0.000382689177
NL: BROYDEN 12 ; 0.0114571067 0.000188258059
NL: BROYDEN 13 ; 0.00554434897 9.11022655e-05
NL: BROYDEN 14 ; 0.00259570587 4.26514792e-05
NL: BROYDEN 15 ; 0.0011364197 1.86731408e-05
NL: BROYDEN 16 ; 0.000434046425 7.13205696e-06
NL: BROYDEN 17 ; 0.000125657613 2.06474977e-06
NL: BROYDEN 18 ; 2.10704958e-05 3.46220975e-07
NL: BROYDEN 19 ; 1.29557101e-06 2.12882441e-08
NL: BROYDEN 20 ; 1.45277217e-08 2.38713033e-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.3887968 0.367882708
NL: BROYDEN 2 ; 12.510447 0.205566076
NL: BROYDEN 3 ; 5.98619238 0.098362439
NL: BROYDEN 4 ; 3.04366868 0.050012204
NL: BROYDEN 5 ; 1.51175582 0.0248404963
NL: BROYDEN 6 ; 0.75759176 0.0124484094
NL: BROYDEN 7 ; 0.378237612 0.00621503147
NL: BROYDEN 8 ; 0.188999522 0.00310555572
NL: BROYDEN 9 ; 0.0942955644 0.00154942259
NL: BROYDEN 10 ; 0.0469601147 0.000771627627
NL: BROYDEN 11 ; 0.0232898966 0.000382689177
NL: BROYDEN 12 ; 0.0114571067 0.000188258059
NL: BROYDEN 13 ; 0.00554434897 9.11022655e-05
NL: BROYDEN 14 ; 0.00259570587 4.26514792e-05
NL: BROYDEN 15 ; 0.0011364197 1.86731408e-05
NL: BROYDEN 16 ; 0.000434046425 7.13205696e-06
NL: BROYDEN 17 ; 0.000125657613 2.06474977e-06
NL: BROYDEN 18 ; 2.10704958e-05 3.46220975e-07
NL: BROYDEN 19 ; 1.29557101e-06 2.12882441e-08
NL: BROYDEN 20 ; 1.45277217e-08 2.38713033e-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.

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.3887968 0.367882708
NL: BROYDEN 2 ; 8.23629948 0.135335194
NL: BROYDEN 3 ; 3.02978539 0.04978408
NL: BROYDEN 4 ; 1.11437824 0.0183109654
NL: BROYDEN 5 ; 0.409725119 0.00673242012
NL: BROYDEN 6 ; 0.150492273 0.0024728218

NL: BROYDEN 7 ; 0.0551239641 0.000905772353
NL: BROYDEN 8 ; 0.0200403261 0.000329293686
NL: BROYDEN 9 ; 0.00713760492 0.000117281936
NL: BROYDEN 10 ; 0.00240281896 3.94820478e-05
NL: BROYDEN 11 ; 0.000692025597 1.13710555e-05
NL: BROYDEN 12 ; 0.000282811849 4.64703796e-06
NL: BROYDEN 13 ; 6.22844815e-05 1.02343078e-06
NL: BROYDEN 14 ; 7.64241042e-06 1.25576674e-07
NL: BROYDEN 15 ; 2.45096367e-07 4.02731402e-09
NL: BROYDEN 16 ; 1.00716012e-09 1.65492052e-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 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()

promotes=['*'])

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.