最小化具有 2 个约束的函数时出现问题 - Python

Question

我正在编写一个程序来最小化受约束和边界约束的几个参数的函数。以防万一你想运行程序，函数由：

def Fnret(mins):
    Bj, Lj, a, b = mins.reshape(4,N)
    S1 = 0; S2 = 0
    Binf = np.zeros(N); Linf = np.zeros(N);   
    for i in range(N):
        sbi=(Bi/2); sli=(Li/2)
        for j in range(i+1):
            sbi -= Bj[j]
            sli -= Lj[j]
        Binf[i]=sbi
        Linf[i]=sli
    
    for i in range(N):
        S1 += (C*(1-np.sin(a[i]))+T*np.sin(a[i])) * ((2*Bj[i]*Binf[i]+Bj[i]**2)/(np.tan(b[i])*np.cos(a[i]))) +\
            (C*(1-np.sin(b[i]))+T*np.sin(b[i])) * ((2*Bj[i]*Linf[i]+Lj[i]*Bj[i])/(np.sin(b[i])))
    S2 += (gamma*Bj[0]/(6*np.tan(b[0])))*((Bi/2)*(Li/2) + 4*(Binf[0]+Bj[0])*(Linf[0]+Lj[0]) + Binf[0]*Linf[0])
    j=1
    while j<(N):
        S2 += (gamma*Bj[j]/(6*np.tan(b[j])))*(Binf[j-1]*Linf[j-1] + 4*(Binf[j]+Bj[j])*(Linf[j]+Lj[j]) + Binf[j]*Linf[j])
        j += 1
    F = 2*(S1+S2)
    return F

其中Bj、Lj、a和b是N维数组给出的最小化结果，N是输入程序，我仔细检查了功能，它工作正常。我的限制条件是：

def Brhs(mins):  # Constraint over Bj
    return np.sum(mins.reshape(4,N)[0]) - Bi

def Lrhs(mins): # Constraint over Lj
    return np.sum(mins.reshape(4,N)[1]) - Li

cons = [{'type': 'eq', 'fun': lambda Bj: 1.0*Brhs(Bj)},
         {'type': 'eq', 'fun': lambda Lj: 1.0*Lrhs(Lj)}]

这样 Bj 的所有分量之和必须等于 Bi（与 Lj 相同）。变量的边界由下式给出：

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

为了使问题可重现，使用以下输入很重要：

# Inputs:
gamma = 17.     
C = 85.        
T = C         
Li = 1.        
Bi = 0.5           
N = 3

为了最小化，我使用了 cyipopt 库（与 scipy.optimize 类似）。然后，最小化由下式给出：

from cyipopt import minimize_ipopt
x0 = np.ones(4*N) # Initial guess 
res = minimize_ipopt(Fnret, x0, constraints=cons, bounds=bounds)

问题是结果不符合我对约束施加的条件（即 Bj 或 Lj 值的总和与结果上的 Bi 或 Li 不同）。但是，例如，如果我只写两个约束之一（在 Lj 或 Bj 上），它对那个变量就可以正常工作。也许我在使用 2 个约束时遗漏了一些东西并且我找不到错误，它似乎不能同时使用两个约束。任何帮助将不胜感激。提前致谢！

P.S.: 此外，我希望函数结果 F 也为正。我怎样才能强加这个条件？谢谢！

Answer 1

不是一个完整的答案，只是一些任意顺序的提示：

您的初始点 x0 不可行，因为它与您的两个约束相矛盾。通过在 x0 评估您的约束，可以很容易地观察到这一点。在幕后，Ipopt 通常会检测到这一点并尝试找到可行的初始点。但是，强烈建议尽可能提供可行的初始点。
您的变量范围不是 well-defined。在您的范围内评估您的 objective 会产生多个除以零的结果。例如，当且仅当 b[i] = 0 时分母 np.tan(b[i]) 为零，因此 0 不是所有 b[i] 的有效值。对其他项进行类似的操作，您应该获得 0 < b[i] < pi/2 和 0 ≤ a[i] < pi/2。在这里，您可以通过 0 + eps ≤ b[i] ≤ pi/2 - eps 和 0 ≤ a[i] ≤ pi/2 - eps 对严格的不等式进行建模，其中 eps 是一个足够小的正数。
如果你真的想强加 objective 函数总是正的，你可以简单地添加不等式约束 Fnret(x) >= 0，即 {'type': 'ineq', 'fun': Fnret}.

在代码中：

# bounds
eps = 1e-8
bounds = [(0, None)]*2*N + [(0, np.pi/2 - eps)]*N + [(0+eps, np.pi/2 - eps)]*N

# (feasible) initial guess
x0 = eps*np.ones(4*N)
x0[[0, N]] = [Bi-(N-1)*eps, Li-(N-1)*eps]

# constraints
cons = [{'type': 'eq', 'fun': Brhs},
        {'type': 'eq', 'fun': Lrhs},
        {'type': 'ineq', 'fun': Fnret}]

res = minimize_ipopt(Fnret, x0, constraints=cons, bounds=bounds, options={'disp': 5})

最后但同样重要的是，这仍然没有收敛到一个固定点，所以很可能确实没有局部最小值。从这里开始，您可以尝试使用其他（可行的！）初始点和 double-check 您问题的数学。还值得提供确切的梯度和约束 Jacobians。

Answer 2

所以，根据@joni 的建议，我可以通过采用 scipy.optimize.minimize 库的 trust-constr 方法找到一个遵守约束的驻点。我的代码是运行如下：

import numpy as np
from scipy.optimize import minimize

# Inputs:
gamma = 17      
C = 85.        
T = C         
Li = 2.         
Bi = 1.           
N = 3          # for instance

# Constraints:
def Brhs(mins):  
    return np.sum(mins.reshape(4,N)[0]) - Bi/2
def Lrhs(mins):
    return np.sum(mins.reshape(4,N)[1]) - Li/2

# Function to minimize:
def Fnret(mins):
    Bj, Lj, a, b = mins.reshape(4,N)
    S1 = 0; S2 = 0
    Binf = np.zeros(N); Linf = np.zeros(N);   
    for i in range(N):
        sbi=(Bi/2); sli=(Li/2)
        for j in range(i+1):
            sbi -= Bj[j]
            sli -= Lj[j]
        Binf[i]=sbi
        Linf[i]=sli
    
    for i in range(N):
        S1 += (C*(1-np.sin(a[i]))+T*np.sin(a[i])) * ((2*Bj[i]*Binf[i]+Bj[i]**2)/(np.tan(b[i])*np.cos(a[i]))) +\
            (C*(1-np.sin(b[i]))+T*np.sin(b[i])) * ((2*Bj[i]*Linf[i]+Lj[i]*Bj[i])/(np.sin(b[i])))
    S2 += (gamma*Bj[0]/(6*np.tan(b[0])))*((Bi/2)*(Li/2) + 4*(Binf[0]+Bj[0])*(Linf[0]+Lj[0]) + Binf[0]*Linf[0])
    j=1
    while j<(N):
        S2 += (gamma*Bj[j]/(6*np.tan(b[j])))*(Binf[j-1]*Linf[j-1] + 4*(Binf[j]+Bj[j])*(Linf[j]+Lj[j]) + Binf[j]*Linf[j])
        j += 1
    F = 2*(S1+S2)
    return F

eps = 1e-3
bounds = [(0,None)]*2*N + [(0+eps,np.pi/2-eps)]*2*N       # Bounds

cons = ({'type': 'ineq', 'fun': Fnret},
        {'type': 'eq', 'fun':  Lrhs},
        {'type': 'eq', 'fun':  Brhs})

x0 = np.ones(4*N) # Initial guess
res = minimize(Fnret, x0, method='trust-constr', bounds = bounds, constraints=cons, tol=1e-6)

F = res.fun
Bj = (res.x).reshape(4,N)[0]
Lj = (res.x).reshape(4,N)[1]
ai = (res.x).reshape(4,N)[2]
bi = (res.x).reshape(4,N)[3]

本质上是一样的，只是改变了最小化技术。从 np.sum(Bj) 和 np.sum(Lj) 很容易看出，结果与施加的约束一致，这在以前是行不通的。

最小化具有 2 个约束的函数时出现问题 - Python

Problems minimizing a function with 2 constraints - Python

python

minimization

scipy

scipy-optimize

cyipopt