从矩阵动态构建scipy.optimize.minimize的目标函数和约束

Question

在我的代码中，我正在尝试写一个循环，这个循环会使用scipy.optimize.minimize中的SQSLP选项来更新x的值，因为我有一些不等式和等式的约束。除了更新x，这个过程还会更新目标函数和约束中的值。

下面是调用所有函数的部分，这些函数会在给定的x值下更新约束和目标函数。

import numpy as np
from scipy.optimize import minimize

num_vars = 2
num_eq_cons = 1
num_ineq_cons = 1


# return the desired initial value for x
def x_init():
    x = np.zeros([num_vars,1], dtype = float)
    x[0] = 2
    x[1] = 1
    return x

# Objective
# needs the current value of x (x_k)
def eval_objective(x):
    f = x[0]**2 + x[1]**2
    return float(f)

# Inequality Constraints
# need the current value of x (x_k)
def eval_ineq_cons(x):
    g = np.zeros([num_ineq_cons,1], dtype = float)
    g[0] = x[0]**2 - x[1]
    return g

# Equality Constraints
# need the current value of x (x_k)
def eval_eq_cons(x):
    h = np.zeros([num_eq_cons, 1], dtype = float)
    h[0] = 2 - x[0] - x[1]**2
    return h

# Taylor series Approximations
# needs the current value of x
# For objective
def eval_part_obj(x):
    dfdx = np.zeros([num_vars,1], dtype = float)
    dfdx[0] = 2*x[0]
    dfdx[1] = 2*x[1]
    return dfdx

# For inequalities
def eval_part_ineq_cons(x):
    dgdx = np.zeros([num_ineq_cons, num_vars], dtype = float)
    dgdx[0,0] = 2*x[0]
    dgdx[0,1] = -1.0
    return dgdx

# For equalities
def eval_part_eq_cons(x):
    dhdx = np.zeros([num_eq_cons, num_vars], dtype = float)
    dhdx[0,0] = -1
    dhdx[0,1] = -2*x[1]
    return dhdx

x_k = x_init() # Initialization point
f_k = eval_objective(x_k)
dfdx = eval_part_obj(x_k)
g_k = eval_ineq_cons(x_k)
dgdx = eval_part_ineq_cons(x_k)
h_k = eval_eq_cons(x_k)
dhdx = eval_part_eq_cons(x_k)

下面的代码是我出错的地方：

fun = lambda x_k1: f_k + dfdx[0]*(x_k1[0] - x_k[0]) + dfdx[1]*(x_k1[1] - x_k[1])

cons = ({'type': 'ineq',
         'fun': lambda x_k1: g_k + np.dot(dgdx, (x_k1 - x_k)).squeeze(),
         'jac': dgdx},
        {'type': 'eq',
         'fun': lambda x_k1: h_k + np.dot(dhdx, (x_k1 - x_k)).squeeze(),
         'jac': dhdx})

bnds = [(.5, 2.5), (0,3)]

res = minimize( fun, x_k, method = 'SLSQP', bounds=bnds, constraints = cons)
print(res)

我得到的错误是：

Traceback (most recent call last):
  File "SLPmain.py", line 104, in <module>
    res = minimize( fun, x_k, method = 'SLSQP', bounds=bnds, constraints = cons)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/_minimize.py", line 358, in minimize
    constraints, **options)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/slsqp.py", line 376, in _minimize_slsqp
    for con in cons['eq']])
TypeError: 'numpy.ndarray' object is not callable

如果有更好的方法来编写约束和目标函数，请告诉我！

Answer 1

这里有两个明显的问题：

1. 不等式和等式约束的雅可比矩阵（Jacobian）都是 array 类型，而不是可以调用的函数。它们的定义是：

dgdx = eval_part_ineq_cons(x_k)
dhdx = eval_part_eq_cons(x_k)

它们是函数的返回值（因此是 array），而不是函数本身。

2. 约束函数的返回值不是标量（单个数值），而是 array。例如：

g_k + np.dot(dgdx, (x_k1 - x_k)).squeeze()

要解决这个问题，具体方法取决于你的实际工作问题，有一些替代方案，比如：

(g_k + np.dot(dgdx, (x_k1 - x_k)).squeeze()).sum()

或者

((g_k + np.dot(dgdx, (x_k1 - x_k)).squeeze())**2).sum() #sum of squares.

不等式约束也有同样的问题。我怀疑这是因为 fun 被作为 eqcons 传递给了 fmin_slsqp，而不是 f_eqcons。如果你希望等式约束函数的所有元素返回值都等于0，可能应该直接使用 fmin_slsqp 并传入约束函数 f_eqcons。

因此，可以这样做：

cons = ({'type': 'ineq',
         'fun': lambda x_k1: (g_k + np.dot(dgdx, (x_k1 - x_k)).squeeze()).sum(),
         'jac':  eval_part_ineq_cons},
        {'type': 'eq',
         'fun': lambda x_k1: (h_k + np.dot(dhdx, (x_k1 - x_k)).squeeze()).sum(),
         'jac': eval_part_eq_cons})

结果是：

In [31]:

res
Out[31]:
  status: 9
 success: False
    njev: 101
    nfev: 1385
     fun: array([ 2.69798824])
       x: array([ 1.89963627,  0.04972158])
 message: 'Iteration limit exceeded'
     jac: array([ 4.,  2.,  0.])
     nit: 101

实际上，注意到优化失败了。我怀疑当你有雅可比约束时，它可能不会保证成功。去掉雅可比约束，它就能工作（虽然可能不是你想要的方式）。

In [29]:

res
Out[29]:
  status: 0
 success: True
    njev: 11
    nfev: 122
     fun: array([ 1.33333333])
       x: array([ 0.5       ,  2.16666667])
 message: 'Optimization terminated successfully.'
     jac: array([ 4.,  2.,  0.])
     nit: 15