Source code for openmdao.lib.drivers.conmindriver

"""
    ``conmindriver.py`` - Driver for the CONMIN optimizer.
    
    The CONMIN driver can be a good example of how to wrap another Driver for
    OpenMDAO. However, there are some things to keep in mind.
    
    1. This implementation of the CONMIN Fortran driver is interruptable, in
    that control is returned every time an objective or constraint evaluation
    is needed. Most external optimizers just expect a function to be passed
    for these, so they require a slightly different driver implementation than
    this.
    
    2. The CONMIN driver is a direct wrap, and all inputs are passed using
    Numpy arrays. This means there are a lot of locally stored variables
    that you might not need for a pure Python optimizer.
    
    Ultimately, if you are wrapping a new optimizer for OpenMDAO, you should
    endeavour to understand the purpose for each statement so that your
    implementation doesn't do any unneccessary or redundant calculation.
"""

# pylint: disable-msg=C0103

#public symbols
__all__ = ['CONMINdriver']

import logging

# pylint: disable-msg=E0611,F0401
try:
    from numpy import zeros, ones
    from numpy import int as numpy_int
    import conmin.conmin as conmin
except ImportError as err:
    logging.warn("In %s: %r" % (__file__, err))
    # to keep class decl from barfing before being stubbed out
    zeros = lambda *args, **kwargs: None 
    numpy_int = int

from openmdao.main.driver_uses_derivatives import DriverUsesDerivatives
from openmdao.main.exceptions import RunStopped
from openmdao.main.datatypes.api import Array, Bool, Enum, Float, Int
from openmdao.main.interfaces import IHasParameters, IHasIneqConstraints, \
                                     IHasObjective, implements, IOptimizer
from openmdao.main.hasparameters import HasParameters
from openmdao.main.hasconstraints import HasIneqConstraints
from openmdao.main.hasobjective import HasObjective
from openmdao.util.decorators import add_delegate, stub_if_missing_deps


class _cnmn1(object):
    """Just a primitive data structure for storing cnmnl common block data.
    
    We save the common blocks to prevent collision in the case where there are
    multiple instances of CONMIN running in our model."""
    
    def __init__(self):
        self.clear()

    def clear(self):
        """ Clear values. """
        
        # pylint: disable-msg=W0201
        self.ndv = 0
        self.ncon = 0
        self.nside = 0
        self.infog = 0
        self.info = 0
        self.nfdg = 0
        self.icndir = 0
        self.nscal = 0
        self.fdch = 0.0
        self.fdchm = 0.0
        self.ct = 0.
        self.ctmin = 0.
        self.ctlmin = 0.
        self.theta = 0.
        self.delfun = 0.
        self.dabfun = 1.e-8
        self.linobj = 0
        self.itrm = 0
        self.alphax = 0.
        self.abobj1 = 0.
        self.ctl = 0.
        self.igoto = 0
        self.nac = 0
        self.iter = 0
        self.iprint = 0
        self.itmax = 0
        # pylint: enable-msg=W0201
 
class _consav(object):
    """Just a primitive data structure for storing consav common block data.
    
    We save the common blocks to prevent collision in the case where there are
    multiple instances of CONMIN running in our model."""
    
    def __init__(self):
        self.clear()
    
    def clear(self):
        """ Clear values. """
        
        # pylint: disable-msg=W0201
        self.dm1 = 0.0
        self.dm2 = 0.0
        self.dm3 = 0.0
        self.dm4 = 0.0
        self.dm5 = 0.0
        self.dm6 = 0.0
        self.dm7 = 0.0
        self.dm8 = 0.0
        self.dm9 = 0.0
        self.dm10 = 0.0
        self.dm11 = 0.0
        self.dm12 = 0.0
        self.dct = 0.0
        self.dctl = 0.0
        self.phi = 0.0
        self.abobj = 0.0
        self.cta = 0.0
        self.ctam = 0.0
        self.ctbm = 0.0
        self.obj1 = 0.0
        self.slope = 0.0
        self.dx = 0.0
        self.dx1 = 0.0
        self.fi = 0.0
        self.xi = 0.0
        self.d_objtdf1 = 0.0
        self.alp = 0.0
        self.fff = 0.0
        self.a1 = 0.0
        self.a2 = 0.0
        self.a3 = 0.0
        self.a4 = 0.0
        self.f1 = 0.0
        self.f2 = 0.0
        self.f3 = 0.0
        self.f4 = 0.0
        self.cv1 = 0.0
        self.cv2 = 0.0
        self.cv3 = 0.0
        self.cv4 = 0.0
        self.app = 0.0
        self.alpca = 0.0
        self.alpfes = 0.0
        self.alpln = 0.0
        self.alpmin = 0.0
        self.alpnc = 0.0
        self.alpsav = 0.0
        self.alpsid = 0.0
        self.alptot = 0.0
        self.rspace = 0.0
        self.idm1 = 0
        self.idm2 = 0
        self.idm3 = 0
        self.jdir = 0
        self.iobj = 0
        self.kobj = 0
        self.kcount = 0
        self.ncal = [0, 0]
        self.nfeas = 0
        self.mscal = 0
        self.ncobj = 0
        self.nvc = 0
        self.kount = 0
        self.icount = 0
        self.igood1 = 0
        self.igood2 = 0
        self.igood3 = 0
        self.igood4 = 0
        self.ibest = 0
        self.iii = 0
        self.nlnc = 0
        self.jgoto = 0
        self.ispace = [0, 0]
        # pylint: enable-msg=W0201

@stub_if_missing_deps('numpy', 'conmin')
@add_delegate(HasParameters, HasIneqConstraints, HasObjective)
[docs]class CONMINdriver(DriverUsesDerivatives): """ Driver wrapper of Fortran version of CONMIN. Note on self.cnmn1.igoto, which reports CONMIN's operation state: 0: Initial and final state 1: Initial evaluation of objective and constraint values 2: Evalute gradients of objective and constraints (internal only) 3: Evalute gradients of objective and constraints 4: One-dimensional search on unconstrained function 5: Solve 1D search problem for unconstrained function """ # I don't see an IUsesGradients implements(IHasParameters, IHasIneqConstraints, IHasObjective, IOptimizer) # pylint: disable-msg=E1101 # Control parameters for CONMIN. # CONMIN has quite a few parameters to give the user control over aspects # of the solution. cons_is_linear = Array(zeros(0,'i'), dtype=numpy_int, iotype='in', desc='Array designating whether each constraint is linear.') iprint = Enum(0, [0, 1, 2, 3, 4, 5, 101], iotype='in', desc='Print ' 'information during CONMIN solution. Higher values are ' 'more verbose. 0 suppresses all output.') itmax = Int(10, iotype='in', desc='Maximum number of iterations before ' 'termination.') fdch = Float(.01, iotype='in', desc='Relative change in parameters ' 'when calculating finite difference gradients.' ' (only when CONMIN calculates gradient)') fdchm = Float(.01, iotype='in', desc='Minimum absolute step in finite ' 'difference gradient calculations.' ' (only when CONMIN calculates gradient)') icndir = Float(0, iotype='in', desc='Conjugate gradient restart. ' 'parameter.') ct = Float(-0.1, iotype='in', desc='Constraint thickness parameter.') ctmin = Float(0.004, iotype='in', desc='Minimum absolute value of ct ' 'used in optimization.') ctl = Float(-0.01, iotype='in', desc='Constraint thickness parameter for ' 'linear and side constraints.') ctlmin = Float(0.001, iotype='in', desc='Minimum absolute value of ctl ' 'used in optimization.') theta = Float(1.0, iotype='in', desc='Mean value of the push-off factor ' 'in the method of feasible directions.') phi = Float(5.0, iotype='in', desc='Participation coefficient - penalty ' 'parameter that pushes designs towards the feasible ' 'region.') delfun = Float(0.001, iotype='in', low=0.0, desc='Relative convergence tolerance.') dabfun = Float(0.001, iotype='in', low=1.0e-10, desc='Absolute convergence tolerance.') linobj = Bool(False, iotype='in', desc='Linear objective function flag. ' 'Set to True if objective is linear.') itrm = Int(3, iotype='in', desc='Number of consecutive iterations to ' 'indicate convergence (relative or absolute).') def __init__(self, *args, **kwargs): super(CONMINdriver, self).__init__(*args, **kwargs) # Save data from common blocks into our CONMINdriver object self.cnmn1 = _cnmn1() self.consav = _consav() self.iter_count = 0 # Since CONMIN is an external Fortran code, it requires data to be # passed via several Numpy arrays. We definie all these arrays here # and initialize to zero length. They are later resized in # config_conmin. # basic stuff self.design_vals = zeros(0,'d') self._scal = zeros(0,'d') self.cons_active_or_violated = zeros(0, 'i') self._cons_is_linear = zeros(0, 'i') self.d_const = zeros(0, 'd') # gradient of objective w.r.t x[i] self.d_obj = zeros(0, 'd') # move direction in the optimization space self.s = zeros(0, 'd') # temp storage self._b = zeros(0, 'd') self._c = zeros(0, 'd') self._ms1 = zeros(0, 'i') # temp storage for constraints self.g1 = zeros(0,'d') self.g2 = zeros(0,'d')
[docs] def start_iteration(self): """Perform initial setup before iteration loop begins.""" # Flag used to figure out if we are starting a new finite difference self.baseline_point = True self._config_conmin() self.cnmn1.igoto = 0 self.iter_count = 0 # get the initial values of the parameters # check if any min/max constraints are violated by initial values for i, val in enumerate(self.get_parameters().values()): self.design_vals[i] = dval = val.evaluate(self.parent) vlow = val.low vhigh = val.high if dval > vhigh: if (dval - vhigh) < self.ctlmin: self.design_vals[i] = vhigh else: self.raise_exception('initial value of: %s is greater than maximum' % val.target, ValueError) if dval < vlow: if (vlow - dval) < self.ctlmin: self.design_vals[i] = vlow else: self.raise_exception('initial value of: %s is less than minimum' % val.target, ValueError)
[docs] def continue_iteration(self): """Returns True if iteration should continue.""" return self.cnmn1.igoto != 0 or self.iter_count == 0
[docs] def pre_iteration(self): """Checks for RunStopped.""" super(CONMINdriver, self).pre_iteration() if self._stop: self.raise_exception('Stop requested', RunStopped)
[docs] def run_iteration(self): """ The CONMIN driver iteration.""" #self._logger.debug('iter_count = %d' % self.iter_count) #self._logger.debug('objective = %f' % self.cnmn1.obj) #self._logger.debug('design vals = %s' % self.design_vals[:-2]) # TODO: 'step around' ill-behaved cases. self._load_common_blocks() #print "Iteration %s: " % self.get_pathname(), self.iter_count #print "Before" #print self.design_vals try: (self.design_vals, self._scal, self.d_const, self.s, self.g1, self.g2, self._b, self._c, self._cons_is_linear, self.cons_active_or_violated, self._ms1) = \ conmin.conmin(self.design_vals, self._lower_bounds, self._upper_bounds, self.constraint_vals, self._scal, self.d_obj, self.d_const, self.s, self.g1, self.g2, self._b, self._c, self._cons_is_linear, self.cons_active_or_violated, self._ms1) except Exception, err: self._logger.error(str(err)) raise self._save_common_blocks() # calculate objective and constraints if self.cnmn1.info == 1: # Note. CONMIN is driving the finite difference estimation of the # gradient. However, we still take advantage of a component's # user-defined gradients via Fake Finite Difference. if self.cnmn1.igoto == 3: # Save baseline states and calculate derivatives if self.baseline_point: self.calc_derivatives(first=True) self.baseline_point = False # update the parameters in the model dvals = [float(val) for val in self.design_vals[:-2]] self.set_parameters(dvals) # Run model under Fake Finite Difference self.ffd_order = 1 super(CONMINdriver, self).run_iteration() self.ffd_order = 0 else: # update the parameters in the model dvals = [float(val) for val in self.design_vals[:-2]] self.set_parameters(dvals) # Run the model for this step super(CONMINdriver, self).run_iteration() self.baseline_point = True # calculate objective self.cnmn1.obj = self.eval_objective() # update constraint value array for i, v in enumerate(self.get_ineq_constraints().values()): val = v.evaluate(self.parent) if '>' in val[2]: self.constraint_vals[i] = val[1]-val[0] else: self.constraint_vals[i] = val[0]-val[1] #self._logger.debug('constraints = %s'%self.constraint_vals) # calculate gradient of constraints and gradient of objective # We also have to determine which constrints are active/violated, and # only return gradients of active/violated constraints. elif self.cnmn1.info == 2 and self.cnmn1.nfdg == 1: self.ffd_order = 1 self.differentiator.calc_gradient() self.ffd_order = 0 self.d_obj[:-2] = self.differentiator.get_gradient(self.get_objectives().keys()[0]) for i in range(len(self.cons_active_or_violated)): self.cons_active_or_violated[i] = 0 self.cnmn1.nac = 0 for i, name in enumerate(self.get_ineq_constraints().keys()): if self.constraint_vals[i] >= self.cnmn1.ct: self.cons_active_or_violated[self.cnmn1.nac] = i+1 self.d_const[:-2, self.cnmn1.nac] = \ self.differentiator.get_gradient(name) self.cnmn1.nac += 1 else: self.raise_exception('Unexpected value for flag INFO returned \ from CONMIN.', RuntimeError)
[docs] def post_iteration(self): """ Checks CONMIN's return status and writes out cases.""" super(CONMINdriver, self).post_iteration() # Iteration count comes from CONMIN. You can't just count over the # loop because some cycles do other things (e.g., numerical # gradient calculation) if (self.iter_count != self.cnmn1.iter) or \ self.cnmn1.igoto == 0: self.iter_count = self.cnmn1.iter self.record_case()
def _config_conmin(self): """Set up arrays for the Fortran conmin routine, perform some validation, and make sure that array sizes are consistent. """ self.cnmn1.clear() self.consav.clear() params = self.get_parameters().values() # size arrays based on number of parameters num_dvs = len(params) self.design_vals = zeros(num_dvs+2, 'd') if num_dvs < 1: self.raise_exception('no parameters specified', RuntimeError) # create lower_bounds array self._lower_bounds = zeros(num_dvs+2) for i, param in enumerate(params): self._lower_bounds[i] = param.low # create upper bounds array self._upper_bounds = zeros(num_dvs+2) for i, param in enumerate(params): self._upper_bounds[i] = param.high # create array for CONMIN's internal scaling # we no longer use these, but the still need the empty arrays self._scal = ones(num_dvs+2) self.s = zeros(num_dvs+2, 'd') # size constraint related arrays length = len(self.get_ineq_constraints()) + 2*num_dvs self.constraint_vals = zeros(length, 'd') # temp storage of constraint and design vals self.g1 = zeros(length, 'd') # temp storage of constraint vals self.g2 = zeros(length, 'd') # if constraint i is known to be a linear function of design vars, # the user can set cons_is_linear[i] to 1, otherwise set it to 0. This # is not essential and is for efficiency only. self._cons_is_linear = zeros(length, 'i') if len(self.cons_is_linear) > 0: if len(self.cons_is_linear) != len(self.get_ineq_constraints()): self.raise_exception('size of cons_is_linear (%d) does not \ match number of constraints (%d)'% (len(self.cons_is_linear),length), ValueError) else: for i, val in enumerate(self.cons_is_linear): self._cons_is_linear[i] = val self.cnmn1.ndv = num_dvs self.cnmn1.ncon = len(self.get_ineq_constraints()) self.cnmn1.nside = 2*num_dvs nacmx1 = max(num_dvs, self.cnmn1.ncon+self.cnmn1.nside)+1 n1 = num_dvs+2 n3 = nacmx1 n4 = max(n3, num_dvs) n5 = 2*n4 # array of active or violated constraints (ic in CONMIN) self.cons_active_or_violated = zeros(n3, 'i') # stuff for gradients if self.differentiator: self.nfdg = 1 else: self.nfdg = 0 self.d_const = zeros((int(n1), int(n3)), 'd') self.d_obj = zeros(num_dvs+2, 'd') # these are all temp storage self._b = zeros((int(n3), int(n3)), 'd') self._c = zeros(n4, 'd') self._ms1 = zeros(n5, 'i') # Load all of the user-changeable parameters into the common block self.cnmn1.iprint = self.iprint self.cnmn1.itmax = self.itmax self.cnmn1.fdch = self.fdch self.cnmn1.fdchm = self.fdchm self.cnmn1.icndir = self.icndir self.cnmn1.nscal = 0 self.cnmn1.nfdg = self.nfdg self.cnmn1.ct = self.ct self.cnmn1.ctmin = self.ctmin self.cnmn1.ctl = self.ctl self.cnmn1.ctlmin = self.ctlmin self.cnmn1.theta = self.theta self.consav.phi = self.phi self.cnmn1.dabfun = self.dabfun self.cnmn1.delfun = self.delfun if self.linobj: self.cnmn1.linobj = 1 else: self.cnmn1.linobj = 0 self.cnmn1.itrm = self.itrm def _load_common_blocks(self): """ Reloads the common blocks using the intermediate info saved in the class. """ for name, value in self.cnmn1.__dict__.items(): setattr( conmin.cnmn1, name, value ) for name, value in self.consav.__dict__.items(): setattr( conmin.consav, name, value ) def _save_common_blocks(self): """" Saves the common block data to the class to prevent trampling by other instances of CONMIN. """ common = self.cnmn1 for name, value in common.__dict__.items(): setattr(common, name, type(value)(getattr(conmin.cnmn1, name))) consav = self.consav for name, value in consav.__dict__.items(): setattr(consav, name, type(value)(getattr(conmin.consav, name)))
OpenMDAO Home