我有一个相当简单的约束优化问题，但会根据我的操作方式得到不同的答案。让我们先完成导入和 pretty-print 功能:

import numpy as np
from scipy.optimize import minimize, LinearConstraint, NonlinearConstraint, SR1

def print_res( res, label ):
    print("\n\n ***** ", label, " ***** \n")
    print(res.message)
    print("obj func value at solution", obj_func(res.x))
    print("starting values: ", x0)
    print("ending values:   ", res.x.astype(int) )
    print("% diff", (100.*(res.x-x0)/x0).astype(int) )
    print("target achieved?",target,res.x.sum())

示例数据非常简单:

n = 5
x0 = np.arange(1,6) * 10_000
target = x0.sum() + 5_000   # increase sum from 15,000 to 20,000

这是约束优化(包括 jacobians)。换句话说，我要最小化的目标函数只是从初始值到最终值的平方百分比变化之和。线性相等 约束只是要求x.sum() 等于一个常数。

def obj_func(x):
    return ( ( ( x - x0 ) / x0 ) ** 2 ).sum()

def obj_jac(x):
    return 2. * ( x - x0 ) / x0 ** 2

def constr_func(x):
    return x.sum() - target

def constr_jac(x):
    return np.ones(n)

为了比较，我通过使用等式约束将 x[0] 替换为 x[1:] 的函数，将其重构为无约束最小化>。请注意，未约束函数传递的是 x0[1:]，而约束函数传递的是 x0。

def unconstr_func(x):
    x_one       = target - x.sum()
    first_term  = ( ( x_one - x0[0] ) / x0[0] ) ** 2
    second_term = ( ( ( x - x0[1:] ) / x0[1:] ) ** 2 ).sum()
    return first_term + second_term

然后我尝试通过三种方式最小化:

1. 不受“Nelder-Mead”的约束
2. 受“trust-constr”约束(w/& w/o jacobian)
3. 受“SLSQP”约束(w/& w/o jacobian)

代码:

##### (1) unconstrained

res0 = minimize( unconstr_func, x0[1:], method='Nelder-Mead')   # OK, but weird note
res0.x = np.hstack( [target - res0.x.sum(), res0.x] )
print_res( res0, 'unconstrained' )    

##### (2a) constrained -- trust-constr w/ jacobian

nonlin_con = NonlinearConstraint( constr_func, 0., 0., constr_jac )
resTCjac = minimize( obj_func, x0, method='trust-constr',
                     jac='2-point', hess=SR1(), constraints = nonlin_con )
print_res( resTCjac, 'trust-const w/ jacobian' )

##### (2b) constrained -- trust-constr w/o jacobian

nonlin_con = NonlinearConstraint( constr_func, 0., 0. )    
resTC = minimize( obj_func, x0, method='trust-constr',
                  jac='2-point', hess=SR1(), constraints = nonlin_con )    
print_res( resTC, 'trust-const w/o jacobian' )

##### (3a) constrained -- SLSQP w/ jacobian

eq_cons = { 'type': 'eq', 'fun' : constr_func, 'jac' : constr_jac }
resSQjac = minimize( obj_func, x0, method='SLSQP',
                     jac = obj_jac, constraints = eq_cons )    
print_res( resSQjac, 'SLSQP w/ jacobian' )

##### (3b) constrained -- SLSQP w/o jacobian

eq_cons = { 'type': 'eq', 'fun' : constr_func }    
resSQ = minimize( obj_func, x0, method='SLSQP',
                  jac = obj_jac, constraints = eq_cons )
print_res( resSQ, 'SLSQP w/o jacobian' )

这是一些简化的输出(当然您可以运行代码以获得完整输出):

starting values:  [10000 20000 30000 40000 50000]

***** (1) unconstrained  *****
Optimization terminated successfully.
obj func value at solution 0.0045454545454545305
ending values:    [10090 20363 30818 41454 52272]

***** (2a) trust-const w/ jacobian  *****
The maximum number of function evaluations is exceeded.
obj func value at solution 0.014635854609684874
ending values:    [10999 21000 31000 41000 51000]

***** (2b) trust-const w/o jacobian  *****
`gtol` termination condition is satisfied.
obj func value at solution 0.0045454545462939935
ending values:    [10090 20363 30818 41454 52272]

***** (3a) SLSQP w/ jacobian  *****
Optimization terminated successfully.
obj func value at solution 0.014636111111111114
ending values:    [11000 21000 31000 41000 51000]    

***** (3b) SLSQP w/o jacobian  *****   
Optimization terminated successfully.
obj func value at solution 0.014636111111111114
ending values:    [11000 21000 31000 41000 51000]

注意事项:

1. (1) 和 (2b) 是合理的解决方案，因为它们实现了显着较低的目标函数值，并且凭直觉我们期望具有较大起始值的变量移动得更多(绝对值和百分比)较小的。

2. 将 jacobian 添加到“trust-const”会导致它得到错误的答案(或至少是更糟糕的答案)并且还会超过最大迭代次数。也许 jacobian 是错误的，但功能非常简单，我很确定它是正确的 (?)

3. 'SLSQP' 似乎在没有或没有提供的 jacobian 的情况下工作，但工作速度非常快，并声称可以成功终止。这似乎非常令人担忧，因为得到错误的答案并声称已成功终止几乎是最糟糕的结果。

4. 最初，我使用了非常小的起始值和目标(仅为上述值的 1/1,000)，在这种情况下，上述所有 5 种方法都可以正常工作并给出相同的答案。我的示例数据仍然非常小，处理 1,2,..,5 而不是 1000,2000,..5000 似乎有点奇怪。

5. FWIW，请注意，3 个不正确的结果都通过向每个初始值增加 1,000 来达到目标​​ - 这满足约束条件，但远未达到最小化目标函数(具有较高初始值的 b/c 变量应该是增加的多于减少的，以最小化百分比差异的平方和)。

所以我的问题实际上就是这里发生了什么，为什么只有 (1) 和 (2b) 似乎有效？

更一般地说，我想找到一个很好的基于 python 的方法来解决这个和类似的优化问题，并且会考虑使用除 scipy 之外的其他包的答案，尽管最好的答案最好也解决 scipy 在这里发生的事情(例如这是用户错误还是我应该发布到 github 的错误？)。

最佳答案

下面是如何使用 nlopt 解决这个问题，它是一个非线性优化库，给我留下了深刻的印象。

首先，目标函数和梯度都使用相同的函数定义:

def obj_func(x, grad):
    if grad.size > 0:
        grad[:] = obj_jac(x)
    return ( ( ( x/x0 - 1 )) ** 2 ).sum()

def obj_jac(x):
    return 2. * ( x - x0 ) / x0 ** 2

def constr_func(x, grad):
    if grad.size > 0:
        grad[:] = constr_jac(x)
    return x.sum() - target

def constr_jac(x):
    return np.ones(n)

然后，使用 Nelder-Mead 和 SLSQP 运行最小化:

opt = nlopt.opt(nlopt.LN_NELDERMEAD,len(x0)-1)
opt.set_min_objective(unconstr_func)
opt.set_ftol_abs(1e-15)
xopt = opt.optimize(x0[1:].copy())
xopt = np.hstack([target - xopt.sum(), xopt])
fval = opt.last_optimum_value()
print_res(xopt,fval,"Nelder-Mead");

opt = nlopt.opt(nlopt.LD_SLSQP,len(x0))
opt.set_min_objective(obj_func)
opt.add_equality_constraint(constr_func)
opt.set_ftol_abs(1e-15)
xopt = opt.optimize(x0.copy())
fval = opt.last_optimum_value()
print_res(xopt,fval,"SLSQP w/ jacobian");

结果如下:

 *****  Nelder-Mead  ***** 

obj func value at solution 0.00454545454546
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10090 20363 30818 41454 52272]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0


 *****  SLSQP w/ jacobian  ***** 

obj func value at solution 0.00454545454545
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10090 20363 30818 41454 52272]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0

在对此进行测试时，我想我发现了最初尝试的问题所在。如果我将函数的绝对公差设置为 1e-8，这是 scipy 函数的默认值，我得到:

 *****  Nelder-Mead  ***** 

obj func value at solution 0.0045454580693
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10090 20363 30816 41454 52274]
% diff [0 1 2 3 4]
target achieved? 155000.0 155000.0


 *****  SLSQP w/ jacobian  ***** 

obj func value at solution 0.0146361108503
result:  3
starting values:  [ 10000.  20000.  30000.  40000.  50000.]
ending values:    [10999 21000 31000 41000 51000]
% diff [9 5 3 2 2]
target achieved? 155000.0 155000.0

这正是您所看到的。所以我的猜测是，在 SLSQP 期间，最小化器最终出现在似然空间中的某个位置，其中下一个跳跃距离最后一个位置小于 1e-8。

关于python - 使用 SciPy 最小化受线性等式约束的二次函数，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/55543140/