I'm using COBYLA to do a cost minimization on a linear objective function with constraints. I'm implementing lower and upper bounds by including a constraint for each.
import numpy as np
import scipy.optimize
def linear_cost(factor_prices):
def cost_fn(x):
return np.dot(factor_prices, x)
return cost_fn
def cobb_douglas(factor_elasticities):
def tech_fn(x):
return np.product(np.power(x, factor_elasticities), axis=1)
return tech_fn
def mincost(targets, cost_fn, tech_fn, bounds):
n = len(bounds)
m = len(targets)
x0 = np.ones(n) # Do not use np.zeros.
cons = []
for factor in range(n):
lower, upper = bounds[factor]
l = {'type': 'ineq',
'fun': lambda x: x[factor] - lower}
u = {'type': 'ineq',
'fun': lambda x: upper - x[factor]}
cons.append(l)
cons.append(u)
for output in range(m):
t = {'type': 'ineq',
'fun': lambda x: tech_fn(x)[output] - targets[output]}
cons.append(t)
res = scipy.optimize.minimize(cost_fn, x0,
constraints=cons,
method='COBYLA')
return res
COBYLA doesn't respect the upper or lower-bound constraints, but it does respect the technology constraint.
>>> p = np.array([5., 20.])
>>> cost_fn = linear_cost(p)
>>> fe = np.array([[0.5, 0.5]])
>>> tech_fn = cobb_douglas(fe)
>>> bounds = [[0.0, 15.0], [0.0, float('inf')]]
>>> mincost(np.array([12.0]), cost_fn, tech_fn, bounds)
x: array([ 24.00010147, 5.99997463])
message: 'Optimization terminated successfully.'
maxcv: 1.9607782064667845e-10
nfev: 75
status: 1
success: True
fun: 239.99999999822359
Why wouldn't COBYLA respect the first factor constraint (i.e. upper-bound @ 15)?
COBYLA is in fact respecting all the bounds you give.
The problem lies in the construction of the cons
list.
Namely, binding of variables in lambda and other inner-scoped functions in Python (and Javascript) is lexical, and does not behave in the way you assume: http://eev.ee/blog/2011/04/24/gotcha-python-scoping-closures/ After the loop is finished, the variables lower
and upper
have values 0
and inf
, and variable factor
has value 1
, and these values are what is then used by all of the lambda functions.
One workaround is to explicitly bind the specific values of variables to dummy keyword arguments:
for factor in range(n):
lower, upper = bounds[factor]
l = {'type': 'ineq',
'fun': lambda x, a=lower, i=factor: x[i] - a}
u = {'type': 'ineq',
'fun': lambda x, b=upper, i=factor: b - x[i]}
cons.append(l)
cons.append(u)
for output in range(m):
t = {'type': 'ineq',
'fun': lambda x, i=output: tech_fn(x)[i] - targets[i]}
cons.append(t)
A second way is to add a factory function generating the lambdas.