使用scipy最小化多变量函数.衍生物未知

Question

我有一个函数实际上是对另一个程序的调用(一些Fortran代码).当我调用这个函数(run_moog)时,我可以解析4个变量,并返回6个值.这些值都应该接近0(为了最小化).但是,我把它们组合起来:np.sum(results**2).现在我有一个标量函数.我想最小化这个功能,即np.sum(results**2)尽可能接近零.
注意:当此函数(run_moog)接受4个输入参数时,它会为Fortran代码创建一个依赖于这些参数的输入文件.

我已经尝试了几种方法来从scipy docs中优化它.但没有一个按预期工作.最小化应该能够对4个变量进行限制.这是一个尝试:

from scipy.optimize import minimize # Tried others as well from the docs
x0 = 4435, 3.54, 0.13, 2.4
bounds = [(4000, 6000), (3.00, 4.50), (-0.1, 0.1), (0.0, None)]
a = minimize(fun_mmog, x0, bounds=bounds, method='L-BFGS-B')  # I've tried several different methods here
print a

然后这给了我

  status: 0
 success: True
    nfev: 5
     fun: 2.3194639999999964
       x: array([  4.43500000e+03,   3.54000000e+00,   1.00000000e-01,
         2.40000000e+00])
 message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     jac: array([ 0., 0., -54090399.99999981, 0.])
     nit: 0

第三个参数略有变化,而其他参数完全相同.还有5个函数调用(nfev)但没有迭代(nit).scipy的输出显示在这里.

Answer 1

几种可能性:

1. 试试COBYLA.它应该是无衍生的,并支持不等式约束.
2. 你不能通过普通的接口使用不同的epsilons; 所以尝试用1e4缩放你的第一个变量.(将它分开,然后再多出来.)
3. 跳过正常的自动jacobian构造函数,并自己创建:

假设您正在尝试使用SLSQP,并且您没有提供雅可比功能.它适合你.它的代码approx_jacobian在slsqp.py中.这是一个浓缩版本:

def approx_jacobian(x,func,epsilon,*args):
    x0 = asfarray(x)
    f0 = atleast_1d(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()

你可以尝试用以下方法替换该循环:

    for (i, e) in zip(range(len(x0)), epsilon):
        dx[i] = e
        jac[i] = (func(*((x0+dx,)+args)) - f0)/e
        dx[i] = 0.0

你不能把它作为jacobian提供minimize,但修复它是很简单的:

def construct_jacobian(func,epsilon):
    def jac(x, *args):
        x0 = asfarray(x)
        f0 = atleast_1d(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()
    return jac

然后你可以这样打电话minimize:

minimize(fun_mmog, x0,
         jac=construct_jacobian(fun_mmog, [1e0, 1e-4, 1e-4, 1e-4]),
         bounds=bounds, method='SLSQP')