Source code for openmdao.lib.drivers.broydensolver

"""
    ``broyednsolver.py`` -- Solver based on the nonlinear solvers found in ``Scipy.Optimize``.
       
"""

# pylint: disable-msg=C0103

#public symbols
__all__ = ['BroydenSolver']

import logging

try:
    import numpy
except ImportError as err:
    logging.warn("In %s: %r" % (__file__, err))
else:
    # this little funct replaces a dependency on scipy
    npnorm = numpy.linalg.norm
    def norm(a, ord=None):
        return npnorm(numpy.asarray_chkfinite(a), ord=ord)

# pylint: disable-msg=E0611,F0401
from openmdao.lib.datatypes.api import Float, Int, Enum
                                 
from openmdao.main.api import Driver
from openmdao.main.exceptions import RunStopped
from openmdao.main.hasparameters import HasParameters
from openmdao.main.hasconstraints import HasEqConstraints
from openmdao.util.decorators import add_delegate, stub_if_missing_deps
from openmdao.main.interfaces import IHasParameters, IHasEqConstraints, \
                                     ISolver, implements

    
@stub_if_missing_deps('numpy')
@add_delegate(HasParameters, HasEqConstraints)
[docs]class BroydenSolver(Driver): """ :term:`MIMO` Newton-Raphson Solver with Broyden approximation to the Jacobian. Algorithms are based on those found in ``scipy.optimize``. Nonlinear solvers These solvers find x for which F(x)=0. Both x and F are multidimensional. All solvers have the parameter *iter* (the number of iterations to compute); some of them have other parameters of the solver. See the particular solver for details. A collection of general-purpose nonlinear multidimensional solvers. - ``broyden2``: Broyden's second method -- the same as ``broyden1`` but updates the inverse Jacobian directly - ``broyden3``: Broyden's second method -- the same as ``broyden2`` but instead of directly computing the inverse Jacobian, it remembers how to construct it using vectors. When computing ``inv(J)*F``, it uses those vectors to compute this product, thus avoiding the expensive NxN matrix multiplication. - ``excitingmixing``: The excitingmixing algorithm. ``J=-1/alpha`` The ``broyden2`` is the best. For large systems, use ``broyden3``; excitingmixing is also very effective. The remaining nonlinear solvers from SciPy are, in their own words, of "mediocre quality," so they were not implemented. """ implements(IHasParameters, IHasEqConstraints, ISolver) # pylint: disable-msg=E1101 algorithm = Enum('broyden2', ['broyden2', 'broyden3', 'excitingmixing'], iotype = 'in', desc='Algorithm to use. Choose from ' 'broyden2, broyden3, and excitingmixing.') itmax = Int(10, iotype='in', desc='Maximum number of iterations before ' 'termination.') alpha = Float(0.4, iotype='in', desc='Mixing Coefficient.') alphamax = Float(1.0, iotype='in', desc='Maximum Mixing Coefficient (only ' 'used with excitingmixing.)') tol = Float(0.00001, iotype='in', desc='Convergence tolerance. If the norm of the independent ' 'vector is lower than this, then terminate successfully.') def __init__(self): super(BroydenSolver, self).__init__() self.xin = numpy.zeros(0,'d') self.F = numpy.zeros(0,'d')
[docs] def execute(self): """Solver execution.""" # get the initial values of the independents independents = self.get_parameters().values() self.xin = numpy.zeros(len(independents),'d') for i, val in enumerate(independents): self.xin[i] = val.evaluate(self.parent) # perform an initial run for self-consistency self.pre_iteration() self.run_iteration() self.post_iteration() # get initial dependents dependents = self.get_eq_constraints().values() self.F = numpy.zeros(len(dependents),'d') for i, val in enumerate(dependents): term = val.evaluate(self.parent) self.F[i] = term[0] - term[1] # pick solver algorithm if self.algorithm == 'broyden2': self.execute_broyden2() elif self.algorithm == 'broyden3': self.execute_broyden3() elif self.algorithm == 'excitingmixing': self.execute_excitingmixing()
[docs] def execute_broyden2(self): """From SciPy, Broyden's second method. Updates inverse Jacobian by an optimal formula. There is NxN matrix multiplication in every iteration. The best norm(F(x))=0.003 achieved in ~20 iterations. """ xm = self.xin.T Fxm = numpy.matrix(self.F).T Gm = -self.alpha*numpy.matrix(numpy.identity(len(self.xin))) for n in range(self.itmax): if self._stop: self.raise_exception('Stop requested', RunStopped) deltaxm = -Gm*Fxm xm = xm + deltaxm.T # update the new independents in the model self.set_parameters(xm.flat) # run the model self.pre_iteration() self.run_iteration() self.post_iteration() # get dependents for i, val in enumerate(self.get_eq_constraints().values()): term = val.evaluate(self.parent) self.F[i] = term[0] - term[1] self.record_case() # successful termination if independents are below tolerance if norm(self.F) < self.tol: return Fxm1 = numpy.matrix(self.F).T deltaFxm = Fxm1 - Fxm if norm(deltaFxm) == 0: msg = "Broyden iteration has stopped converging. Change in " + \ "input has produced no change in output. This could " + \ "indicate a problem with your component connections. " + \ "It could also mean that this solver method is " + \ "inadequate for your problem." raise RuntimeError(msg) Fxm = Fxm1.copy() Gm = Gm + (deltaxm-Gm*deltaFxm)*deltaFxm.T/norm(deltaFxm)**2
[docs] def execute_broyden3(self): """from scipy, Broyden's second (sic) method. Updates inverse Jacobian by an optimal formula. The NxN matrix multiplication is avoided. The best norm(F(x))=0.003 achieved in ~20 iterations. """ zy = [] def updateG(z, y): """G:=G+z*y.T'""" zy.append((z, y)) def Gmul(f): """G=-alpha*1+z*y.T+z*y.T ...""" s = -self.alpha*f for z, y in zy: s = s + z*(y.T*f) return s xm = self.xin.T Fxm = numpy.matrix(self.F).T for n in range(self.itmax): if self._stop: self.raise_exception('Stop requested', RunStopped) deltaxm = Gmul(-Fxm) xm = xm + deltaxm.T # update the new independents in the model self.set_parameters(xm.flat) # run the model self.pre_iteration() self.run_iteration() self.post_iteration() # get dependents for i, val in enumerate(self.get_eq_constraints().values()): term = val.evaluate(self.parent) self.F[i] = term[0] - term[1] self.record_case() # successful termination if independents are below tolerance if norm(self.F) < self.tol: return Fxm1 = numpy.matrix(self.F).T deltaFxm = Fxm1 - Fxm if norm(deltaFxm) == 0: msg = "Broyden iteration has stopped converging. Change in " + \ "input has produced no change in output. This could " + \ "indicate a problem with your component connections. " + \ "It could also mean that this solver method is " + \ "inadequate for your problem." raise RuntimeError(msg) Fxm = Fxm1.copy() updateG(deltaxm - Gmul(deltaFxm), deltaFxm/norm(deltaFxm)**2)
[docs] def execute_excitingmixing(self): """from scipy, The excitingmixing method. J=-1/alpha The best norm(F(x))=0.005 achieved in ~140 iterations. Note: SciPy uses 0.1 as the default value for alpha for this algorithm. Ours is set at 0.4, which is appropriate for ``Broyden2`` and ``Broyden3``, so adjust it accordingly if there are problems. """ xm = self.xin.copy() beta = numpy.array([self.alpha]*len(xm)) Fxm = self.F.T.copy() for n in range(self.itmax): if self._stop: self.raise_exception('Stop requested', RunStopped) deltaxm = beta*Fxm xm = xm + deltaxm # update the new independents in the model self.set_parameters(xm.flat) # run the model self.pre_iteration() self.run_iteration() self.post_iteration() # get dependents for i, val in enumerate(self.get_eq_constraints().values()): term = val.evaluate(self.parent) self.F[i] = term[0] - term[1] self.record_case() # successful termination if independents are below tolerance if norm(self.F) < self.tol: return Fxm1 = self.F.T for i in range(len(xm)): if Fxm1[i]*Fxm[i] > 0: beta[i] = beta[i] + self.alpha if beta[i] > self.alphamax: beta[i] = self.alphamax else: beta[i] = self.alpha Fxm = Fxm1.copy() # end broydensolver.py
OpenMDAO Home