[Scipy-svn] r3696 - trunk/scipy/optimize
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Dec 22 12:20:14 EST 2007
Author: rob.falck
Date: 2007-12-22 11:20:01 -0600 (Sat, 22 Dec 2007)
New Revision: 3696
Modified:
trunk/scipy/optimize/slsqp.py
Log:
slsqp updated to allow constraints via functions fixed and additional bounds checking
Modified: trunk/scipy/optimize/slsqp.py
===================================================================
--- trunk/scipy/optimize/slsqp.py 2007-12-21 02:04:21 UTC (rev 3695)
+++ trunk/scipy/optimize/slsqp.py 2007-12-22 17:20:01 UTC (rev 3696)
@@ -1,19 +1,55 @@
-__all__ = ['fmin_slsqp']
+"""This module implements the Sequential Least SQuares Programming optimization
+algorithm (SLSQP), orginally developed by Dieter Kraft.
+See http://www.netlib.org/toms/733
+
+"""
+
+__all__ = ['approx_jacobian','fmin_slsqp']
+
from _slsqp import slsqp
from numpy import zeros, array, identity, linalg, rank, squeeze, append, \
- asfarray,product
-from sys import exit
-from math import sqrt
+ asfarray,product, concatenate, finfo, sqrt, vstack, transpose
from optimize import approx_fprime, wrap_function
-import numpy
-_epsilon = sqrt(numpy.finfo(float).eps)
+__docformat__ = "restructuredtext en"
+_epsilon = sqrt(finfo(float).eps)
+
+def approx_jacobian(x,func,epsilon,*args):
+ """Approximate the Jacobian matrix of callable function func
+
+ *Parameters*:
+ x - The state vector at which the Jacobian matrix is desired
+ func - A vector-valued function of the form f(x,*args)
+ epsilon - The peturbation used to determine the partial derivatives
+ *args - Additional arguments passed to func
+
+ *Returns*:
+ An array of dimensions (lenf, lenx) where lenf is the length
+ of the outputs of func, and lenx is the number of
+
+ *Notes*:
+ The approximation is done using forward differences
+
+ """
+ x0 = asfarray(x)
+ f0 = func(*((x0,)+args))
+ jac = zeros([len(x0),len(f0)])
+ dx = zeros(len(x0))
+ for i in range(len(x0)):
+ dx[i] = epsilon
+ jac[i] = (func(*((x0+dx,)+args)) - f0)/epsilon
+ dx[i] = 0.0
+ return jac.transpose()
+
+
+
+
def fmin_slsqp( func, x0 , eqcons=[], f_eqcons=None, ieqcons=[], f_ieqcons=None,
- bounds = [], fprime = None, fprime_cons=None,args = (),
- iter = 100, acc = 1.0E-6, iprint = 1, full_output = 0,
- epsilon = _epsilon ):
+ bounds = [], fprime = None, fprime_eqcons=None,
+ fprime_ieqcons=None, args = (), iter = 100, acc = 1.0E-6,
+ iprint = 1, full_output = 0, epsilon = _epsilon ):
"""
Minimize a function using Sequential Least SQuares Programming
@@ -46,16 +82,16 @@
for each independent variable [(xl0, xu0),(xl1, xu1),...]
fprime : callable f(x,*args)
A function that evaluates the partial derivatives of func.
- fprime_cons : callable f(x,*args)
+ fprime_eqcons : callable f(x,*args)
A function of the form f(x, *args) that returns the m by n
- array of constraint normals. If not provided, the normals
- will be approximated. Equality constraint normals precede
- inequality constraint normals. The array returned by
- fprime_cons should be sized as ( len(eqcons) +
- len(ieqcons), len(x0) ). If fprime_cons is not supplied
- (normals are approximated) then the constraints must be
- supplied via the eqcons and ieqcons structures, not
- f_eqcons and f_ieqcons.
+ array of equality constraint normals. If not provided,
+ the normals will be approximated. The array returned by
+ fprime_eqcons should be sized as ( len(eqcons), len(x0) ).
+ fprime_ieqcons : callable f(x,*args)
+ A function of the form f(x, *args) that returns the m by n
+ array of inequality constraint normals. If not provided,
+ the normals will be approximated. The array returned by
+ fprime_ieqcons should be sized as ( len(ieqcons), len(x0) ).
args : sequence
Additional arguments passed to func and fprime.
iter : int
@@ -114,56 +150,79 @@
7 : "Rank-deficient equality constraint subproblem HFTI",
8 : "Positive directional derivative for linesearch",
9 : "Iteration limit exceeded" }
+
+ # Now do a lot of function wrapping
- # Wrap the functions
# Wrap func
feval, func = wrap_function(func, args)
+ # Wrap fprime, if provided, or approx_fprime if not
if fprime:
- # Wrap fprime, if provided
geval, fprime = wrap_function(fprime,args)
else:
- # Wrap approx_fprime, if fprime not provided
geval, fprime = wrap_function(approx_fprime,(func,epsilon))
- if fprime_cons:
- approx_constraint_norms = False
- if f_eqcons:
- ceval, f_eqcons = wrap_function(f_eqcons,args)
+
+ if f_eqcons:
+ # Equality constraints provided via f_eqcons
+ ceval, f_eqcons = wrap_function(f_eqcons,args)
+ if fprime_eqcons:
+ # Wrap fprime_eqcons
+ geval, fprime_eqcons = wrap_function(fprime_eqcons,args)
else:
- for i in range(len(eqcons)):
- if eqcons[i]:
- ceval, eqcons[i] = wrap_function(eqcons[i],args)
- if f_ieqcons:
- ceval, f_ieqcons = wrap_function(f_ieqcons,args)
- else:
- for i in range(len(ieqcons)):
- if ieqcons[i]:
- ceval, ieqcons[i] = wrap_function(ieqcons[i],args)
- geval, fprime_cons = wrap_function(fprime_cons,args)
+ # Wrap approx_jacobian
+ geval, fprime_eqcons = wrap_function(approx_jacobian,
+ (f_eqcons,epsilon))
else:
- approx_constraint_norms = True
+ # Equality constraints provided via eqcons[]
eqcons_prime = []
for i in range(len(eqcons)):
eqcons_prime.append(None)
if eqcons[i]:
+ # Wrap eqcons and eqcons_prime
ceval, eqcons[i] = wrap_function(eqcons[i],args)
- geval, eqcons_prime[i] = \
- wrap_function(approx_fprime, (eqcons[i],epsilon))
+ geval, eqcons_prime[i] = wrap_function(approx_fprime,
+ (eqcons[i],epsilon))
+
+ if f_ieqcons:
+ # Inequality constraints provided via f_ieqcons
+ ceval, f_ieqcons = wrap_function(f_ieqcons,args)
+ if fprime_ieqcons:
+ # Wrap fprime_ieqcons
+ geval, fprime_ieqcons = wrap_function(fprime_ieqcons,args)
+ else:
+ # Wrap approx_jacobian
+ geval, fprime_ieqcons = wrap_function(approx_jacobian,
+ (f_ieqcons,epsilon))
+ else:
+ # Inequality constraints provided via ieqcons[]
ieqcons_prime = []
for i in range(len(ieqcons)):
ieqcons_prime.append(None)
if ieqcons[i]:
+ # Wrap ieqcons and ieqcons_prime
ceval, ieqcons[i] = wrap_function(ieqcons[i],args)
- geval, ieqcons_prime[i] = \
- wrap_function(approx_fprime, (ieqcons[i],epsilon))
+ geval, ieqcons_prime[i] = wrap_function(approx_fprime,
+ (ieqcons[i],epsilon))
+
# Transform x0 into an array.
x = asfarray(x0).flatten()
# Set the parameters that SLSQP will need
- meq = len(eqcons) # meq = The number of equality constraints
- m = meq + len(ieqcons) # m = The total number of constraints
- la = array([1,m]).max() # la =
- n = len(x) # n = The number of independent variables
+ # meq = The number of equality constraints
+ if f_eqcons:
+ meq = len(f_eqcons(x))
+ else:
+ meq = len(eqcons)
+ if f_ieqcons:
+ mieq = len(f_ieqcons(x))
+ else:
+ mieq = len(ieqcons)
+ # m = The total number of constraints
+ m = meq + mieq
+ # la = The number of constraints, or 1 if there are no constraints
+ la = array([1,m]).max()
+ # n = The number of independent variables
+ n = len(x)
# Define the workspaces for SLSQP
n1 = n+1
@@ -173,14 +232,23 @@
len_jw = mineq
w = zeros(len_w)
jw = zeros(len_jw)
-
+
# Decompose bounds into xl and xu
if len(bounds) == 0:
bounds = [(-1.0E12, 1.0E12) for i in range(n)]
- if len(bounds) != n:
- raise IndexError, 'SLSQP Error: If bounds is specified, len(bounds) == len(x0)'
+ elif len(bounds) != n:
+ raise IndexError, \
+ 'SLSQP Error: If bounds is specified, len(bounds) == len(x0)'
+ else:
+ for i in range(len(bounds)):
+ if bounds[i][0] > bounds[i][1]:
+ raise ValueError, \
+ 'SLSQP Error: lb > ub in bounds[' + str(i) +'] ' + str(bounds[4])
+
xl = array( [ b[0] for b in bounds ] )
- xu = array( [ b[1] for b in bounds ] )
+ xu = array( [ b[1] for b in bounds ] )
+
+
# Initialize the iteration counter and the mode value
mode = array(0,int)
@@ -193,36 +261,73 @@
print "%5s %5s %16s %16s" % ("NIT","FC","OBJFUN","GNORM")
while 1:
+
if mode == 0 or mode == 1: # objective and constraint evaluation requird
+
# Compute objective function
fx = func(x)
# Compute the constraints
if f_eqcons:
- ceq = f_eqcons(x)
+ c_eq = f_eqcons(x)
else:
- ceq = [ eqcons[i](x) for i in range(meq) ]
+ c_eq = array([ eqcons[i](x) for i in range(meq) ])
if f_ieqcons:
- cieq = f_ieqcons(x)
+ c_ieq = f_ieqcons(x)
else:
- cieq = [ ieqcons[i](x) for i in range(len(ieqcons)) ]
- c = numpy.concatenate( (ceq,cieq), 1)
- #c = array ( [ eqcons[i](x) for i in range(meq) ] +
- # [ ieqcons[i](x) for i in range(len(ieqcons)) ] )
+ c_ieq = array([ ieqcons[i](x) for i in range(len(ieqcons)) ])
+
+ # Now combine c_eq and c_ieq into a single matrix
+ if m == 0:
+ # no constraints
+ c = zeros([la])
+ else:
+ # constraints exist
+ if meq > 0 and mieq == 0:
+ # only equality constraints
+ c = c_eq
+ if meq == 0 and mieq > 0:
+ # only inequality constraints
+ c = c_ieq
+ if meq > 0 and mieq > 0:
+ # both equality and inequality constraints exist
+ c = append(c_eq, c_ieq)
+
if mode == 0 or mode == -1: # gradient evaluation required
+
# Compute the derivatives of the objective function
# For some reason SLSQP wants g dimensioned to n+1
g = append(fprime(x),0.0)
- # Compute the normals of the constraints
- if approx_constraint_norms:
- a = zeros([la,n1])
+ # Compute the normals of the constraints
+ if fprime_eqcons:
+ a_eq = fprime_eqcons(x)
+ else:
+ a_eq = zeros([meq,n])
for i in range(meq):
- a[i] = append(eqcons_prime[i](x),0.0)
- for i in range(meq+1,m):
- a[i] = append(ieqcons_prime[i](x),0.0)
+ a_eq[i] = eqcons_prime[i](x)
+
+ if fprime_ieqcons:
+ a_ieq = fprime_ieqcons(x)
else:
- a = numpy.concatenate( (fprime_cons(x),zeros([la,1])),1)
-
+ a_ieq = zeros([mieq,n])
+ for i in range(mieq):
+ a_ieq[i] = ieqcons_prime[i](x)
+
+ # Now combine a_eq and a_ieq into a single a matrix
+ if m == 0:
+ # no constraints
+ a = zeros([la,n])
+ elif meq > 0 and mieq == 0:
+ # only equality constraints
+ a = a_eq
+ elif meq == 0 and mieq > 0:
+ # only inequality constraints
+ a = a_ieq
+ elif meq > 0 and mieq > 0:
+ # both equality and inequality constraints exist
+ a = vstack((a_eq,a_ieq))
+ a = concatenate((a,zeros([la,1])),1)
+
# Call SLSQP
slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw)
@@ -246,7 +351,7 @@
print " Function evaluations:", feval[0]
print " Gradient evaluations:", geval[0]
- if full_output == 0:
+ if not full_output:
return x
else:
return [list(x),
More information about the Scipy-svn
mailing list