使用 scipy.optimize.minimize 最小化 Python 中的约束函数时出现问题

Question

我正在尝试最小化采用算法 scipy.optimize.minimize 的几个变量的约束函数。该函数涉及 3*N 参数的最小化，其中 N 是输入。更具体地说，我的最小化参数在三个数组 H = H[0],H[1],...,H[N-1]、a = a[0],a[1],...,a[N-1] 和 b = b[0],b[1],...,b[N-1] 中给出，我将它们连接在一个名为 mins 和 len(mins)=3*N 的数组中。 =34=]

这些参数还受到如下约束：

0 <= H  and  sum(H) = 0.5
0 <= a <= Pi/2
0 <= b <= Pi/2

所以，我的约束代码如下：

import numpy as np

# constraints on x:
def Hlhs(mins):    # left hand side
    return np.diag(np.ones(N)) @ mins.reshape(3,N)[0]
def Hrhs(mins):    # right hand side
    return np.sum(mins.reshape(3,N)[0]) - 0.5
con1H = {'type': 'ineq', 'fun': lambda H: Hlhs(H)}
con2H = {'type': 'eq', 'fun': lambda H: Hrhs(H)}

# constraints on a:
def alhs(mins):
    return np.diag(np.ones(N)) @ mins.reshape(3,N)[1]
def arhs(mins): 
    return -np.diag(np.ones(N)) @ mins.reshape(3,N)[1] + (np.ones(N))*np.pi/2
con1a = {'type': 'ineq', 'fun': lambda a: alhs(a)}
con2a = {'type': 'ineq', 'fun': lambda a: arhs(a)}

# constraints on b:
def blhs(mins):
    return np.diag(np.ones(N)) @ mins.reshape(3,N)[2]
def brhs(mins): 
    return -np.diag(np.ones(N)) @ mins.reshape(3,N)[2] + (np.ones(N))*np.pi/2
con1b = {'type': 'ineq', 'fun': lambda b: blhs(b)}
con2b = {'type': 'ineq', 'fun': lambda b: brhs(b)}

我的函数，在其他参数（并采用N=3）被最小化后，给出了（如果太长请见谅）：

gamma = 17      
C = 85         
T = 0         
Hf = 0.5         
Li = 2         
Bi = 1  
         
N = 3            

def FUN(mins):
    H, a, b = mins.reshape(3,N)
    S1 = 0; S2 = 0
    B = np.zeros(N); L = np.zeros(N);   
    for i in range(N):
        sbi=Bi; sli=Li
        for j in range(i+1):
            sbi += 2*H[j]*np.tan(b[j])
            sli += 2*H[j]*np.tan(a[j])
        B[i]=sbi
        L[i]=sli    
    for i in range(N):
        S1 += (C*(1-np.sin(a[i])) + T*np.sin(a[i])) * (Bi*H[i]+H[i]**2*np.tan(b[i]))/np.cos(a[i]) + \
        (C*(1-np.sin(b[i])) + T*np.sin(b[i])) * (Li*H[i]+H[i]**2*np.tan(a[i]))/np.cos(b[i])    
    S2 += (gamma*H[0]/12)*(Bi*Li + 4*(B[0]-H[0]*np.tan(b[0]))*(L[0]-H[0]*np.tan(a[0])) + B[0]*L[0])
    j=1 
    while j<(N):
        S2 += (gamma*H[j]/12)*(B[j-1]*L[j-1] + 4*(B[j]-H[j]*np.tan(b[j]))*(L[j]-H[j]*np.tan(a[j])) + B[j]*L[j])
        j += 1
    F = 2*(S1+S2)
    return F

最后，采用初始猜测值 0，最小化由下式给出：

x0 = np.zeros(3*N)
res = scipy.optimize.minimize(FUN,x0,constraints=(con1H,con2H,con1a,con2a,con1b,con2b),tol=1e-25)

我的问题是：

a) 观察结果 res，有些值变成负值，尽管我有限制它们为正值。最小化成功False，消息为：Positive directional derivative for linesearch。此外，结果与预期的最低要求相去甚远。

b) 采用 method='trust-constr' 我得到的值更接近我的预期，但错误的成功和消息 The maximum number of function evaluations is exceeded.。有什么办法可以改善吗？

我知道有一个非常接近这些值的最小值：

H = [0.2,0.15,0.15]
a = [1.0053,1.0053,1.2566]
b = [1.0681,1.1310,1.3195]

其中函数的值为 123,45。我已经多次检查该功能，它似乎工作正常。谁能帮我找到我的问题出在哪里？我尝试更改 xtol 和 maxiter 但没有成功。

Answer 1

这里有一些提示：

你的初始点x0不可行，因为它不满足约束sum(H) = 0.5。提供一个可行的初始点应该可以解决您的第一个问题。
除了约束sum(H) = 0.5，所有约束都是变量的简单边界。一般情况下，建议通过minimize的bounds参数传递变量边界。您可以像这样简单地定义并传递所有边界

from scipy.optimize import minimize
import numpy as np

# ..your variables and functions ..

bounds = [(0, None)]*N + [(0, np.pi/2)]*2*N

x0 = np.zeros(3*N)
x0[0] = 0.5

res = minimize(FUN, x0, constraints=(con2H,), bounds=bounds,
               method="trust-constr", options={'maxiter': 20000})

其中每个元组包含每个变量的下限和上限。

不幸的是，'trust-constr' 仍然难以收敛到局部最小值。在这种情况下，您可以尝试其他初始点，也可以改用 state-of-the-art 开源求解器 Ipopt。 Cython 包装器 cyipopt 提供类似于 scipy:

from cyipopt import minimize_ipopt

# rest as above

res = minimize_ipopt(FUN, x0, constraints=(con2H,), bounds=bounds)

这给了我一个 objective 值为 122.9 的解决方案。

最后但同样重要的是，提供精确的梯度、雅可比矩阵和粗麻布矩阵总是一个好主意。

使用 scipy.optimize.minimize 最小化 Python 中的约束函数时出现问题

Problem minimizing a constrained function in Python with scipy.optimize.minimize

python

optimization

numeric

scipy

scipy-optimize