"""
``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
