在Python中定义scipy优化的动态约束

Question

我想抽象出以下函数，该函数计算目标函数的最小值以及当我们可以获得任意数量的 g 的最小值时的值。

我从两个变量的简单情况开始，效果很好

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: s[0] * x[0] + s[1] * x[1]

    cons = [
        {'type': 'ineq', 'fun': lambda x: x[0] - ((dist[0] + dist[1]) / eff - g_0)},                   # g_1 > (dist[0] + dist[1]) / eff - g_0
        {'type': 'ineq', 'fun': lambda x: x[1] - ((dist[0] + dist[1] + dist[2]) / eff - x[0] - g_0)},  # g_2 > (dist[0] + dist[1] + dist[2]) / eff - g_1 - g_0
        {'type': 'ineq', 'fun': lambda x: g_max - (dist[0] / eff - g_0) - x[0]},                       # g_1 < g_max - (dist[0] / eff - g_0)
        {'type': 'ineq', 'fun': lambda x: g_max - (g_0 + x[0] - (dist[0] - dist[1]) / eff) - x[1]},    # g_2 < g_max - (g_0 + g_1 - (dist[0] - dist[1]) / eff)
    ]

    # General constraints for all g
    for i in range(len(s)):
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    # Bounds for the variables (g_1 and g_2)
    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]

    solution = minimize(objective, x0, method='SLSQP', bounds=[(g1_lower_bound, g1_upper_bound), (0, g_max)], constraints=cons)

    g_1, g_2 = map(round, solution.x)
    return g_1, g_2, round(solution.fun, 2)


g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g1, optimal_g2, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g_1 = {optimal_g1}, g_2 = {optimal_g2}")
print(f"Minimum value of the objective function: {minimum_value}")

然后我开始将其抽象为任意数量的 g (i>=1)。这是不平等的一般规则

到目前为止，我已经看到了这段代码，但是当我尝试发送完全相同的参数时，它给出了不同的结果

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: sum(s[i] * x[i] for i in range(len(x)))

    cons = []
    for i in range(len(s)):
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i] - (sum(dist[:i+2]) / eff - sum([g_0] + x[:i]))})  # g_i > sum of dists from dist[0] to dist[i] / eff - sum of g from g_0 to g_(i-1)

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - (sum([g_0] + x[:i]) - (sum(dist[:i+1]) / eff)) - x[i]})  # g_i < g_max - (sum of g from g_0 to g_(i-1) - sum of dists from dist[0] to dist[i] / eff)
        
        # General constraints to ensure each g is between 0 and g_max
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(objective, x0, method='SLSQP', bounds=[(0, g_max) for _ in range(len(s))], constraints=cons)

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g_values, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g = {optimal_g_values}")
print(f"Minimum value of the objective function: {minimum_value}")

第一个函数给出以下结果，这也是正确的：

最佳值：g_1 = 100，g_2 = 120 目标函数最小值：810.0

但是第二个，在尝试使其成为 g 的任意 len 后，它返回：

最优值：g = [135, 85] 目标函数最小值：862.5

我检查了第二个函数中的 lambda 函数，它们似乎是正确的。另外，当我尝试执行此代码时：

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: sum(s[i] * x[i] for i in range(len(x)))

    cons = [
        {'type': 'ineq', 'fun': lambda x: x[0] - ((dist[0] + dist[1]) / eff - g_0)},                   # g_1 > (dist[0] + dist[1]) / eff - g_0
        {'type': 'ineq', 'fun': lambda x: x[1] - ((dist[0] + dist[1] + dist[2]) / eff - x[0] - g_0)},  # g_2 > (dist[0] + dist[1] + dist[2]) / eff - g_1 - g_0
    ]
    
    for i in range(len(s)):
        #cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i] - (sum(dist[:i+2]) / eff - sum([g_0] + x[:i]))})  # g_i > sum of dists from dist[0] to dist[i] / eff - sum of g from g_0 to g_(i-1)

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - (sum([g_0] + x[:i]) - (sum(dist[:i+1]) / eff)) - x[i]})  # g_i < g_max - (sum of g from g_0 to g_(i-1) - sum of dists from dist[0] to dist[i] / eff)
        
        # General constraints to ensure each g is between 0 and g_max
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(objective, x0, method='SLSQP', bounds=[(0, g_max) for _ in range(len(s))], constraints=cons)

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g_values, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g = {optimal_g_values}")
print(f"Minimum value of the objective function: {minimum_value}")

它给出了正确答案：

最佳值：g_1 = 100，g_2 = 120 目标函数最小值：810.0

但是如果我们申请其他约束：

import numpy as np
from scipy.optimize import minimize

def optimize(g_0, s, g_max, eff, dist):
    objective = lambda x: sum(s[i] * x[i] for i in range(len(x)))

    cons = [
        {'type': 'ineq', 'fun': lambda x: g_max - (dist[0] / eff - g_0) - x[0]},                       # g_1 < g_max - (dist[0] / eff - g_0)
        {'type': 'ineq', 'fun': lambda x: g_max - (g_0 + x[0] - (dist[0] - dist[1]) / eff) - x[1]},    # g_2 < g_max - (g_0 + g_1 - (dist[0] - dist[1]) / eff)
    ]
    
    for i in range(len(s)):
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i] - (sum(dist[:i+2]) / eff - sum([g_0] + x[:i]))})  # g_i > sum of dists from dist[0] to dist[i] / eff - sum of g from g_0 to g_(i-1)

        #cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - (sum([g_0] + x[:i]) - (sum(dist[:i+1]) / eff)) - x[i]})  # g_i < g_max - (sum of g from g_0 to g_(i-1) - sum of dists from dist[0] to dist[i] / eff)
        
        # General constraints to ensure each g is between 0 and g_max
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(objective, x0, method='SLSQP', bounds=[(0, g_max) for _ in range(len(s))], constraints=cons)

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

g_0 = 80
s = [4.5, 3]
g_max = 135
eff = 5
dist = [400, 500, 600]
optimal_g_values, minimum_value = optimize(g_0, s, g_max, eff, dist)
print(f"Optimal values: g = {optimal_g_values}")
print(f"Minimum value of the objective function: {minimum_value}")

它给出了正确的值，但目标函数的最小值是错误的：

最优值：g = [100, 120] 目标函数最小值：810.65

此时我最好的猜测是，要么我在第二个函数中写了错误的 lambda 函数，要么循环中发生了意外的事情

Answer 1

我发现两个错误。一是您最初指定二变量约束的方式存在错误。第二个是N变量约束指定方式的错误。

我还看到了通过使用约束的线性属性来对函数进行总体改进的一些机会。

让我们从原始规范中的错误开始：

        {'type': 'ineq', 'fun': lambda x: g_max - (g_0 + x[0] - (dist[0] - dist[1]) / eff) - x[1]},    # g_2 < g_max - (g_0 + g_1 - (dist[0] - dist[1]) / eff)

根据您尝试实现的不等式，表达式

(dist[0] - dist[1])

应该是

(dist[0] + dist[1])

。

其次，我发现一个与 NumPy 实现方式相关的错误

+

。

如果您使用

+

，对于 Python 列表来说，这意味着“连接”，或者将第二个列表放在第一个列表的末尾。

例如：

>>> [1] + [2, 3]
[1, 2, 3]
>>> [1] + []
[1]
>>> [1, 2] + [3, 4]
[1, 2, 3, 4]

但是，在 NumPy 中，

+

表示添加。如果

+

的任意一个操作数是 NumPy 数组，则这两个数组将相加。

例如：

>>> np.array([1]) + [2, 3]
array([3, 4])
>>> np.array([1]) + []
array([], dtype=float64)
>>> np.array([1, 2]) + [3, 4]
array([4, 6])

这是非常不同的结果。其效果是表达式

sum([g_0] + x[:i])

根据 x 是列表还是数组执行不同的操作。当 SciPy 优化你的函数时，x 将始终是一个 NumPy 数组。

因此，我建议您将

sum([g_0] + x[:i])

替换为

(g_0 + sum(x[:i]))

，这样可以为列表和数组提供相同的结果。

此外，这些约束是多余的：

        cons.append({'type': 'ineq', 'fun': lambda x, i=i: x[i]})                   # g_i > 0
        cons.append({'type': 'ineq', 'fun': lambda x, i=i: g_max - x[i]})           # g_i < g_max

这些与你的边界做同样的事情，并且边界比约束更有效。我会删除这些。

此外，考虑到所有约束都是 x 和一些常量的线性函数，您可能会发现 scipy.optimize.LinearConstraint 更合适。您的约束（冗余约束除外）相当于以下 LinearConstraint：

    A = np.zeros((len(s), len(s)))
    cons_lb = np.zeros(len(s))
    cons_ub = np.zeros(len(s))
    for i in range(len(s)):
        A[i, :i + 1] = 1
        cons_lb[i] = sum(dist[:i+2]) / eff - g_0
        cons_ub[i] = g_max - g_0 + (sum(dist[:i+1]) / eff)
    cons = LinearConstraint(A, cons_lb, cons_ub)

这有很多好处。 （SciPy 不需要使用 LinearConstraint 进行数值微分，并且计算矩阵乘法比计算多个 Python 函数快得多。）

说到微分，您还可以通过提供雅可比并使用 NumPy 计算目标函数来加快速度。对于 2 个变量，这并不重要，但对于 N 个变量，数值微分需要对目标函数进行 N 次调用，因此对于较大的 N 最好避免使用它。

这是改进后的最终代码：

def optimize(g_0, s, g_max, eff, dist):
    s = np.array(s)
    objective = lambda x: np.sum(s * x)
    jac = lambda x: s

    A = np.zeros((len(s), len(s)))
    cons_lb = np.zeros(len(s))
    cons_ub = np.zeros(len(s))
    for i in range(len(s)):
        A[i, :i + 1] = 1
        cons_lb[i] = sum(dist[:i+2]) / eff - g_0
        cons_ub[i] = g_max - g_0 + (sum(dist[:i+1]) / eff)
    cons = LinearConstraint(A, cons_lb, cons_ub)

    g1_lower_bound = max(0, (dist[0] + dist[1]) / eff - g_0)
    g1_upper_bound = min(g_max, g_max - (dist[0] / eff - g_0))

    # Initial guess for the variables
    x0 = [g1_lower_bound, max(0, ((dist[0] + dist[1] + dist[2]) / eff - g1_lower_bound - g_0) + 1)]
    solution = minimize(
        objective,
        x0,
        jac=jac,
        method='SLSQP',
        bounds=[(0, g_max) for _ in range(len(s))],
        constraints=cons
    )

    g_values = list(map(round, solution.x))
    return g_values, round(solution.fun, 2)

这大约快了 3 倍，并修复了表达式

sum([g_0] + x[:i])

的错误。