Source code for openmdao.solvers.nonlinear.broyden

"""
Define the BroydenSolver class.

Based on implementation in Scipy via OpenMDAO 0.8x with improvements based on NPSS solver.
"""
from __future__ import print_function
import warnings
from six.moves import range

import numpy as np

from openmdao.recorders.recording_iteration_stack import Recording
from openmdao.solvers.solver import NonlinearSolver
from openmdao.utils.class_util import overrides_method
from openmdao.utils.general_utils import simple_warning

CITATION = """@ARTICLE{
              Broyden1965ACo,
              AUTHOR = "C. Broyden",
              TITLE = "A Class of Methods for Solving Nonlinear Simultaneous Equations",
              JOURNAL = "Mathematics of Computation",
              VOLUME = "19",
              YEAR = "1965",
              PAGES = "577--593",
              REFERRED = "[Coleman1996SaE]."
              }"""


[docs]class BroydenSolver(NonlinearSolver): """ Broyden solver. Attributes ---------- delta_fxm : ndarray Most recent change in residual vector. delta_xm : ndarray Most recent change in state vector. fxm : ndarray Most recent residual. Gm : ndarray Most recent Jacobian matrix. linear_solver : LinearSolver Linear solver to use for calculating inverse Jacobian. linesearch : NonlinearSolver Line search algorithm. Default is None for no line search. n : int Total length of the states being solved. xm : ndarray Most recent state. _idx : dict Cache of vector indices for each state name. _computed_jacobians : int Number of computed jacobians. _converge_failures : int Number of consecutive iterations that failed to converge to the tol definied in options. _full_inverse : bool When True, Broyden considers the whole vector rather than a list of states. _recompute_jacobian : bool Flag that becomes True when Broyden detects it needs to recompute the inverse Jacobian. """ SOLVER = 'BROYDEN'
[docs] def __init__(self, **kwargs): """ Initialize all attributes. Parameters ---------- **kwargs : dict options dictionary. """ super(BroydenSolver, self).__init__(**kwargs) # Slot for linear solver self.linear_solver = None # Slot for linesearch self.linesearch = None self.cite = CITATION self.n = 0 self._idx = {} self._recompute_jacobian = True self.Gm = None self.xm = None self.fxm = None self.delta_xm = None self.delta_fxm = None self._converge_failures = 0 self._computed_jacobians = 0 # This gets set to True if the user doesn't declare any states. self._full_inverse = False
def _declare_options(self): """ Declare options before kwargs are processed in the init method. """ super(BroydenSolver, self)._declare_options() self.options.declare('alpha', default=0.4, desc="Value to scale the starting Jacobian, which is Identity. This " "option does nothing if you compute the initial Jacobian " "instead.") self.options.declare('compute_jacobian', types=bool, default=True, desc="When True, compute an initial Jacobian, otherwise start " "with Identity scaled by alpha. Further Jacobians may also be " "computed depending on the other options.") self.options.declare('converge_limit', default=1.0, desc="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.") self.options.declare('cs_reconverge', types=bool, default=True, desc='When True, when this driver solves under a complex step, nudge ' 'the Solution vector by a small amount so that it reconverges.') self.options.declare('diverge_limit', default=2.0, desc="Ratio of current residual to previous residual above which the " "Jacobian will be immediately regenerated.") self.options.declare('max_converge_failures', default=3, desc="The number of convergence failures before regenerating the " "Jacobian.") self.options.declare('max_jacobians', default=10, desc="Maximum number of jacobians to compute.") self.options.declare('state_vars', [], desc="List of the state-variable/residuals that " "are to be solved here.") self.options.declare('update_broyden', default=True, desc="Flag controls whether to perform Broyden update to the " "Jacobian. There are some applications where it may be useful " "to turn this off.") self.supports['gradients'] = True self.supports['implicit_components'] = True def _setup_solvers(self, system, depth): """ Assign system instance, set depth, and optionally perform setup. Parameters ---------- system : <System> Pointer to the owning system. depth : int Depth of the current system (already incremented). """ super(BroydenSolver, self)._setup_solvers(system, depth) self._recompute_jacobian = True self._computed_jacobians = 0 self._disallow_discrete_outputs() if self.linear_solver is not None: self.linear_solver._setup_solvers(self._system, self._depth + 1) else: self.linear_solver = system.linear_solver if self.linesearch is not None: self.linesearch._setup_solvers(self._system, self._depth + 1) self.linesearch._do_subsolve = True states = self.options['state_vars'] prom = system._var_allprocs_prom2abs_list['output'] # Check names of states. bad_names = [name for name in states if name not in prom] if len(bad_names) > 0: msg = "{}: The following variable names were not found: {}" raise ValueError(msg.format(self.msginfo, ', '.join(bad_names))) # Size linear system outputs = system._outputs if len(states) > 0: n = 0 for name in states: size = len(outputs[name]) self._idx[name] = (n, n + size) n += size else: self._full_inverse = True n = len(outputs._data) self.n = n self.Gm = np.empty((n, n)) self.xm = np.empty((n, )) self.fxm = np.empty((n, )) self.delta_xm = None self.delta_fxm = None if self._full_inverse: # Can only use DirectSolver here. from openmdao.solvers.linear.direct import DirectSolver if not isinstance(self.linear_solver, DirectSolver): msg = "{}: Linear solver must be DirectSolver when solving the full model." raise ValueError(msg.format(self.msginfo, ', '.join(bad_names))) return # Always look for states that aren't being solved so we can warn the user. def sys_recurse(system, all_states): subs = system._subsystems_myproc if len(subs) == 0: # Skip implicit components that appear to solve themselves. from openmdao.core.implicitcomponent import ImplicitComponent if overrides_method('solve_nonlinear', system, ImplicitComponent): return all_states.extend(system._list_states()) else: for subsys in subs: sub_nl = subsys.nonlinear_solver if sub_nl and sub_nl.supports['implicit_components']: continue sys_recurse(subsys, all_states) all_states = [] sys_recurse(system, all_states) all_states = [system._var_abs2prom['output'][name] for name in all_states] missing = set(all_states).difference(states) if len(missing) > 0: msg = "The following states are not covered by a solver, and may have been " + \ "omitted from the BroydenSolver 'state_vars': " msg += ', '.join(sorted(missing)) simple_warning(msg) def _assembled_jac_solver_iter(self): """ Return a generator of linear solvers using assembled jacs. """ if self.linear_solver is not None: for s in self.linear_solver._assembled_jac_solver_iter(): yield s def _set_solver_print(self, level=2, type_='all'): """ Control printing for solvers and subsolvers in the model. Parameters ---------- level : int iprint level. Set to 2 to print residuals each iteration; set to 1 to print just the iteration totals; set to 0 to disable all printing except for failures, and set to -1 to disable all printing including failures. type_ : str Type of solver to set: 'LN' for linear, 'NL' for nonlinear, or 'all' for all. """ super(BroydenSolver, self)._set_solver_print(level=level, type_=type_) if self.linear_solver is not None and type_ != 'NL': self.linear_solver._set_solver_print(level=level, type_=type_) if self.linesearch is not None: self.linesearch._set_solver_print(level=level, type_=type_) def _linearize(self): """ Perform any required linearization operations such as matrix factorization. """ if self.linear_solver is not None: self.linear_solver._linearize() if self.linesearch is not None: self.linesearch._linearize() def _iter_initialize(self): """ Perform any necessary pre-processing operations. Returns ------- float Initial relative error in the user-specified residuals. float Initial absolute error in the user-specified residuals. """ system = self._system if self.options['debug_print']: self._err_cache['inputs'] = self._system._inputs._copy_views() self._err_cache['outputs'] = self._system._outputs._copy_views() # Convert local storage if we are under complex step. if system.under_complex_step: if np.iscomplex(self.xm[0]): self.Gm = self.Gm.astype(np.complex) self.xm = self.xm.astype(np.complex) self.fxm = self.fxm.astype(np.complex) elif np.iscomplex(self.xm[0]): self.Gm = self.Gm.real self.xm = self.xm.real self.fxm = self.fxm.real self._converge_failures = 0 self._computed_jacobians = 0 # Execute guess_nonlinear if specified. system._guess_nonlinear() # When under a complex step from higher in the hierarchy, sometimes the step is too small # to trigger reconvergence, so nudge the outputs slightly so that we always get at least # one iteration of Broyden. if system.under_complex_step and self.options['cs_reconverge']: system._outputs._data += np.linalg.norm(self._system._outputs._data) * 1e-10 # Start with initial states. self.xm = self.get_states() with Recording('Broyden', 0, self): self._solver_info.append_solver() # should call the subsystems solve before computing the first residual self._gs_iter() self._solver_info.pop() self._run_apply() norm = self._iter_get_norm() norm0 = norm if norm != 0.0 else 1.0 return norm0, norm def _iter_get_norm(self): """ Return the norm of only the residuals requested in options. Returns ------- float norm. """ # Need to cache the initial residuals, which is done in this function. fxm = self.get_residuals() if not self._full_inverse: # Use full model residual for driving the main loop convergence. fxm = self._system._residuals._data if self._system.under_complex_step: return (np.sum(fxm.real**2) + np.sum(fxm.imag**2))**0.5 else: return np.sum(fxm**2) ** 0.5 def _single_iteration(self): """ Perform the operations in the iteration loop. """ system = self._system Gm = self._update_inverse_jacobian() fxm = self.fxm delta_xm = -Gm.dot(fxm) if self.linesearch: self._solver_info.append_subsolver() self.set_linear_vector(delta_xm) self.linesearch.solve() xm = self.get_states() self._solver_info.pop() else: # Update the new states in the model. xm = self.xm + delta_xm self.set_states(xm) # Run the model. with Recording('Broyden', 0, self): self._solver_info.append_solver() self._gs_iter() self._solver_info.pop() self._run_apply() fxm1 = fxm.copy() fxm = self.get_residuals() delta_fxm = fxm - fxm1 # States may have been further converged hierarchically. xm = self.get_states() delta_xm = xm - self.xm # Determine whether to update Jacobian. self._recompute_jacobian = False opt = self.options if self._computed_jacobians <= opt['max_jacobians']: converge_ratio = np.linalg.norm(fxm) / np.linalg.norm(fxm1) if converge_ratio > opt['diverge_limit']: self._recompute_jacobian = True elif converge_ratio > opt['converge_limit']: self._converge_failures += 1 if self._converge_failures >= opt['max_converge_failures']: self._recompute_jacobian = True else: self._converge_failures = 0 # Cache for next iteration. self.delta_xm = delta_xm self.delta_fxm = delta_fxm self.fxm = fxm self.xm = xm self.Gm = Gm def _update_inverse_jacobian(self): """ Update the inverse Jacobian for a new Broyden iteration. Returns ------- ndarray Updated inverse Jacobian. """ Gm = self.Gm # Apply the Broyden Update approximation to the previous value of the inverse jacobian. if self.options['update_broyden'] and not self._recompute_jacobian: dfxm = self.delta_fxm fact = np.linalg.norm(dfxm) # Sometimes you can get stuck, particularly when enforcing bounds in a linesearch. Make # sure we don't update in this case because of divide by zero. if fact > self.options['atol']: Gm += np.outer((self.delta_xm - Gm.dot(dfxm)), dfxm * (1.0 / fact**2)) # Solve for total derivatives of user-requested residuals wrt states. elif self.options['compute_jacobian']: if self._full_inverse: Gm = self._compute_full_inverse_jacobian() else: Gm = self._compute_inverse_jacobian() self._computed_jacobians += 1 # Set inverse Jacobian to identity scaled by alpha. # This is the default starting point used by scipy and the general broyden algorithm. else: Gm = np.diag(np.full(self.n, -self.options['alpha'], dtype=Gm.dtype)) return Gm
[docs] def get_states(self): """ Return a vector containing the values of the states specified in options. This is used to get the initial state guesses. Returns ------- ndarray Array containing values of states. """ outputs = self._system._outputs if self._full_inverse: xm = outputs._data.copy() else: states = self.options['state_vars'] xm = self.xm.copy() for name in states: i, j = self._idx[name] xm[i:j] = outputs[name] return xm
[docs] def set_states(self, new_val): """ Set new values for states specified in options. Parameters ---------- new_val : ndarray New values for states. """ outputs = self._system._outputs if self._full_inverse: outputs._data[:] = new_val else: states = self.options['state_vars'] for name in states: i, j = self._idx[name] outputs[name] = new_val[i:j]
[docs] def get_residuals(self): """ Return a vector containing the values of the residuals specified in options. Returns ------- ndarray Array containing values of residuals. """ residuals = self._system._residuals if self._full_inverse: fxm = self.fxm = residuals._data.copy() else: states = self.options['state_vars'] fxm = self.fxm for name in states: i, j = self._idx[name] fxm[i:j] = residuals[name] return fxm
[docs] def set_linear_vector(self, dx): """ Copy values from step into the linear vector for backtracking. Parameters ---------- dx : ndarray Full step in the states for this iteration. """ linear = self._system._vectors['output']['linear'] if self._full_inverse: linear._data[:] = dx else: linear.set_const(0.0) states = self.options['state_vars'] for name in states: i, j = self._idx[name] linear[name] = dx[i:j]
def _compute_inverse_jacobian(self): """ Compute inverse Jacobian for system by doing a linear solve for each state. Returns ------- ndarray New inverse Jacobian. """ # TODO : Consider promoting this capability out into OpenMDAO so other solvers can use the # same code. # TODO : Can do each state in parallel if procs are available. system = self._system states = self.options['state_vars'] d_res = system._vectors['residual']['linear'] d_out = system._vectors['output']['linear'] inv_jac = self.Gm d_res.set_const(0.0) # Disable local fd approx_status = system._owns_approx_jac system._owns_approx_jac = False # Linearize model. ln_solver = self.linear_solver do_sub_ln = ln_solver._linearize_children() my_asm_jac = ln_solver._assembled_jac system._linearize(my_asm_jac, sub_do_ln=do_sub_ln) if my_asm_jac is not None and system.linear_solver._assembled_jac is not my_asm_jac: my_asm_jac._update(system) self._linearize() for wrt_name in states: i_wrt, j_wrt = self._idx[wrt_name] d_wrt = d_res[wrt_name] for j in range(j_wrt - i_wrt): # Increment each variable. d_wrt[j] = 1.0 # Solve for total derivatives. ln_solver.solve(['linear'], 'fwd') # Extract results. for of_name in states: i_of, j_of = self._idx[of_name] inv_jac[i_of:j_of, i_wrt + j] = d_out[of_name] d_wrt[j] = 0.0 # Enable local fd system._owns_approx_jac = approx_status return inv_jac def _compute_full_inverse_jacobian(self): """ Compute inverse Jacobian for entire system vector. Only the DirectSolver is supported here. Returns ------- ndarray New inverse Jacobian. """ system = self._system # Disable local fd approx_status = system._owns_approx_jac system._owns_approx_jac = False # Linearize model. ln_solver = self.linear_solver do_sub_ln = ln_solver._linearize_children() my_asm_jac = ln_solver._assembled_jac system._linearize(my_asm_jac, sub_do_ln=do_sub_ln) if my_asm_jac is not None and system.linear_solver._assembled_jac is not my_asm_jac: my_asm_jac._update(system) inv_jac = self.linear_solver._inverse() # Enable local fd system._owns_approx_jac = approx_status return inv_jac def _mpi_print_header(self): """ Print header text before solving. """ if self.options['iprint'] > 0 and self._system.comm.rank == 0: pathname = self._system.pathname if pathname: nchar = len(pathname) prefix = self._solver_info.prefix header = prefix + "\n" header += prefix + nchar * "=" + "\n" header += prefix + pathname + "\n" header += prefix + nchar * "=" print(header)
[docs] def cleanup(self): """ Clean up resources prior to exit. """ super(BroydenSolver, self).cleanup() if self.linear_solver: self.linear_solver.cleanup() if self.linesearch: self.linesearch.cleanup()