与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, and b< /code> 是 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 上），那么它对于该变量来说效果很好。也许我在使用两个约束时遗漏了一些东西，并且找不到错误，似乎它不适用于两个约束。任何帮助将不胜感激。先感谢您！

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

Answer 1

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

• 您的初始点x0是不可行的，因为它与您的两个约束相矛盾。可以通过在x0上评估约束来轻松观察这一点。在引擎盖下，IPOPT通常会检测到这一点，并试图找到可行的初始点。但是，强烈建议尽可能提供可行的初始点。
• 您的变量范围不是明确的。在界限上评估目标可以使多个部门乘以零。例如，分母np.tan（b [i]）在且仅当b [i] = 0时为零，因此0不是所有b的有效值[i] s。同样，对于其他条款，您应该获得0＆lt; b [i]＆lt; pi/2和0≤a[i]＆lt; pi/2。在这里，您可以通过0 +EPS≤b[I]≤pi/2 -EPS和0≤a[i]≤pi/2 -EPS -EPS ，其中eps是一个足够小的正数。
• 如果您真的想强加目标函数始终是正面的，则可以简单地添加不等式约束fnret（x）＆gt; = 0，即{'type'：'ineq'，'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})

最后但并非最不重要的一点是，这仍然不会收敛到固定点，因此确实没有本地最低限度。从这里开始，您可以尝试使用其他（可行的！）初始点进行实验，并仔细检查问题的数学。还值得提供确切的梯度和约束雅各布人。

Answer 2

因此，基于@joni的建议，我可以通过采用trust> cons scipy.optimize.minimize.minimize library的trust cons方法来找到一个固定点。我的代码正在运行如下：

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）很容易看到结果与施加的约束一致，而这些约束以前不起作用。