scipy 约束最优化问题进入受限区域

我正在尝试使用python和scipy来解决多变量优化问题。让我定义一下我正在工作的环境：

搜索的参数：

和问题本身：

(In我的case logL 函数是复杂的，所以我将用一个普通的函数代替它，产生类似的问题。因此，在这个例子中，我没有完全使用函数参数，但是为了问题的一致性，我包括了那些参数。
我使用以下约定在单个平面数组中存储参数：

这是脚本，应该解决我的问题。

import numpy as np
from scipy import optimize as opt
from pprint import pprint
from typing import List

_d = 2
_tmax = 500.0
_T = [[1,2,3,4,5], [6,7,8,9]]

def logL(args: List[float], T : List[List[float]], tmax : float):
    # simplified - normaly using T in computation, here only to determine dimension
    d = len(T)
    # trivially forcing args to go 'out-of constrains'
    return -sum([(args[2 * i] + args[2 * i + 1] * tmax)**2 for i in range(d)]) 

def gradientForIthDimension(i, d, t_max):
    g = np.zeros(2 * d + 2 * d**2)
    g[2 * i] = 1.0
    g[2 * i + 1] = t_max + 1.0

    return g

def zerosWithOneOnJth(j, l):
    r = [0.0 for _ in range(l)]
    r[j] = 1.0

    return r

new_lin_const = {
    'type': 'ineq',
    'fun' : lambda x: np.array(
        [x[2 * i] + x[2 * i + 1] * (_tmax + 1.0) for i in range(_d)] 
        + [x[j] for j in range(2*_d + 2*_d**2) if j not in [2 * i + 1 for i in range(_d)]]
        ),
    'jac' : lambda x: np.array(
        [gradientForIthDimension(i, _d, _tmax) for i in range(_d)]
         + [zerosWithOneOnJth(j, 2*_d + 2*_d**2) for j in range(2*_d + 2*_d**2) if j not in [2 * i + 1 for i in range(_d)]]
        )
}

最后是优化

logArgs = [2 for _ in range(2 * (_d**2) + 2 * _d)]

# addditional bounds, not mentioned in a problem, but suppose a'priori knowledge

bds = [(0.0, 10.0) for _ in range(2 * (_d**2) + 2 * _d)]
for i in range(_d):
    bds[2*i + 1] = (-10.0, 10.0)

res = opt.minimize(lambda x, args: -logL(x, args[0], args[1]),
constraints=new_lin_const, x0 = logArgs, args=([_T, _tmax]), method='SLSQP', options={'disp': True}, bounds=bds)

但在检查结果时，我得到：

pprint(res)

# fun: 2.2124712864600578e-05

# jac: array([0.00665204, 3.32973738, 0.00665204, 3.32973738, 0.        ,

# 0.        , 0.        , 0.        , 0.        , 0.        ,

# 0.        , 0.        ])

# message: 'Optimization terminated successfully'

# nfev: 40

# nit: 3

# njev: 3

# status: 0

# success: True

# x: array([ 1.66633206, -0.00332601,  1.66633206, -0.00332601,  2.        ,

# 2.        ,  2.        ,  2.        ,  2.        ,  2.        ,

# 2.        ,  2.        ])

特别是：

print(res.x[0] + res.x[1]*(501.0))

# -3.2529534621517087e-13

所以结果超出了约束区域...我试图遵循文档，但对我来说它不起作用。我很乐意听到任何关于什么是错误的建议。

首先，请不要再多次发布同一个问题。这个问题和你的另一个问题基本上是一样的here。下次只要编辑你的问题，而不是发布新的问题。
也就是说，你的代码是不必要的复杂，因为你的优化问题是相当简单的。你的目标应该是阅读你的代码就像阅读数学优化问题一样简单。一个更受欢迎的副作用是，它更容易调试你的代码，然后万一它没有按预期工作。
为此，强烈建议您熟悉numpy及其矢量化操作（正如你在上一个问题的注解中提到的）。例如，你不需要循环来实现你的目标、约束函数或雅可比矩阵。将所有优化变量打包到一个大向量x中是正确的方法。然而，你可以简单地将x再次解包成它的lambda，gamma，alpha和beta组件。这使得你更容易编写函数，也更容易阅读。
好吧，你可以在下面找到一个简化的、有效的实现，而不是我在你的代码中一路切下去。通过计算函数并将输出与代码片段中计算出的函数进行比较，你应该知道你这边出了什么问题。

编辑：看起来大多数scipy.minimize下的算法都无法收敛到局部最小值，同时保持约束的严格可行性。如果你愿意使用其他软件包，我建议你使用最先进的NLP求解器Ipopt。你可以通过cyipopt软件包使用它，感谢它的minimize_ipopt方法，您可以像scipy.optimize.minimize一样使用它：

import numpy as np

# from scipy.optimize import minimize

from cyipopt import minimize_ipopt as minimize

d = 2
tmax = 500.0
N = 2*d + 2*d**2

def logL(x, d, tmax):
    lambda_, gamma, alpha, beta = np.split(x, np.cumsum([d, d, d**2]))
    return np.sum((lambda_ + tmax*gamma)**2)

def con_fun(x, d, tmax):
    # split the packed variable x = (lambda_, gamma, alpha, beta)
    lambda_, gamma, alpha, beta = np.split(x, np.cumsum([d, d, d**2]))
    return lambda_ + (tmax + 1.0) * gamma 

def con_jac(x, d, tmax):
    jac = np.block([np.eye(d), (tmax + 1.0)*np.eye(d), np.zeros((d, 2*d**2))])
    return jac

constr = {
    'type': 'ineq', 
    'fun': lambda x: con_fun(x, d, tmax), 
    'jac': lambda x: con_jac(x, d, tmax)
}

bounds = [(0, 10.0)]*N + [(-10.0, 10.0)]*N + [(0.0, 10.0)]*2*d**2
x0 = np.full(N, 2.0)

res = minimize(lambda x: logL(x, d, tmax), x0=x0, constraints=constr, 
    method='SLSQP', options={'disp': True}, bounds=bounds)

print(res)

收益率


******************************************************************************

This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt

******************************************************************************

     fun: 0.00014085582293562834
    info: {'x': array([ 2.0037865 ,  2.0037865 , -0.00399079, -0.00399079,  2.00700641,
        2.00700641,  2.00700641,  2.00700641,  2.00700641,  2.00700641,
        2.00700641,  2.00700641]), 'g': array([0.00440135, 0.00440135]), 'obj_val': 0.00014085582293562834, 'mult_g': array([-0.01675576, -0.01675576]), 'mult_x_L': array([5.00053270e-08, 5.00053270e-08, 1.00240003e-08, 1.00240003e-08,
       4.99251018e-08, 4.99251018e-08, 4.99251018e-08, 4.99251018e-08,
       4.99251018e-08, 4.99251018e-08, 4.99251018e-08, 4.99251018e-08]), 'mult_x_U': array([1.25309309e-08, 1.25309309e-08, 1.00160027e-08, 1.00160027e-08,
       1.25359789e-08, 1.25359789e-08, 1.25359789e-08, 1.25359789e-08,
       1.25359789e-08, 1.25359789e-08, 1.25359789e-08, 1.25359789e-08]), 'status': 0, 'status_msg': b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'}
 message: b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'
    nfev: 15
     nit: 14
    njev: 16
  status: 0
 success: True
       x: array([ 2.0037865 ,  2.0037865 , -0.00399079, -0.00399079,  2.00700641,
        2.00700641,  2.00700641,  2.00700641,  2.00700641,  2.00700641,
        2.00700641,  2.00700641])

并且在所找到的解处评估约束函数产生

In [17]: print(constr['fun'](res.x))
[0.00440135 0.00440135]

因此，满足约束条件。

scipy 约束最优化问题进入受限区域

1条答案

相关问题

热门标签

最新问答