将 scipy.optimize.minimize 限制为整数值

Posted 2023-03-05

【中文标题】将 scipy.optimize.minimize 限制为整数值

【英文标题】：Restrict scipy.optimize.minimize to integer values

【发布时间】：2017-01-07 07:00:09

【问题描述】：

我正在使用scipy.optimize.minimize 来优化答案只能是整数的现实问题。我当前的代码如下所示：

from scipy.optimize import minimize

def f(x):
    return (481.79/(5+x[0]))+(412.04/(4+x[1]))+(365.54/(3+x[2]))+(375.88/(3+x[3]))+(379.75/(3+x[4]))+(632.92/(5+x[5]))+(127.89/(1+x[6]))+(835.71/(6+x[7]))+(200.21/(1+x[8]))

def con(x):
    return sum(x)-7

cons = 'type':'eq', 'fun': con

print scipy.optimize.minimize(f, [1,1,1,1,1,1,1,0,0], constraints=cons, bounds=([0,7],[0,7],[0,7],[0,7],[0,7],[0,7],[0,7],[0,7],[0,7]))

这会产生：

x: array([  2.91950510e-16,   2.44504019e-01,   9.97850733e-01,
     1.05398840e+00,   1.07481251e+00,   2.60570253e-01,
     1.36470363e+00,   4.48527831e-02,   1.95871767e+00]

但我希望它使用整数值进行优化（将所有 x 舍入到最接近的整数并不总是给出最小值）。

有没有办法将scipy.optimize.minimize 仅用于整数值？

（我想我可以创建一个包含x 的所有可能排列的数组，并为每个组合评估 f(x)，但这似乎不是一个非常优雅或快速的解决方案。）

【问题讨论】：

这是不可能的。 numpy/scipy 中没有 (Mixed-)Integer-Programming 的求解器。您可能想要使用pulp 或一些替代品（pyomo、cvxpy、...）。或者如果你疯了：编写你自己的分支定界过程。

您不能将x 舍入为f 内的整数值，并将优化后的结果舍入为整数吗？

【参考方案1】：

纸浆溶液

经过一些研究，我认为您的目标函数不是线性的。我在 Python pulp 库中重新创建了该问题，但纸浆不喜欢我们除以浮点数和“LpAffineExpression”。 This answer 建议线性规划“不理解除法”，但该评论是在添加约束的上下文中，而不是目标函数。该评论将我指向“Mixed Integer Linear Fractional Programming (MILFP)”和Wikipedia。

如果它真的有效，你可以在纸浆中做到这一点（也许有人可以找出原因）：

import pulp

data = [(481.79, 5), (412.04, 4), (365.54, 3)] #, (375.88, 3), (379.75, 3), (632.92, 5), (127.89, 1), (835.71, 6), (200.21, 1)]
x = pulp.LpVariable.dicts('x', range(len(data)), lowBound=0, upBound=7, cat=pulp.LpInteger)

numerator = dict((i,tup[0]) for i,tup in enumerate(data))
denom_int = dict((i,tup[1]) for i,tup in enumerate(data))

problem = pulp.LpProblem('Mixed Integer Linear Programming', sense=pulp.LpMinimize)

# objective function (doesn't work)
# TypeError: unsupported operand type(s) for /: 'float' and 'LpAffineExpression'
problem += sum([numerator[i] / (denom_int[i] + x[i]) for i in range(len(data))])

problem.solve()

for v in problem.variables():
  print(v.name, "=", v.varValue)

---

使用 scipy.optimize 的暴力解决方案

您可以在函数中为每个x 使用brute 和slices 范围。如果您的函数中有 3 个 xs，那么您的范围元组中也将有 3 个 slices。所有这一切的关键是将1 的步长 大小添加到slice(start, stop,step) 所以slice(#, #, 1)。

from scipy.optimize import brute
import itertools

def f(x):
  return (481.79/(5+x[0]))+(412.04/(4+x[1]))+(365.54/(3+x[2]))

ranges = (slice(0, 9, 1),) * 3
result = brute(f, ranges, disp=True, finish=None)
print(result)

---

itertools 解决方案

或者你可以使用 itertools 生成所有组合：

combinations = list(itertools.product(*[[0,1,2,3,4,5,6,7,8]]*3))

values = []
for combination in combinations:
  values.append((combination, f(combination)))

best = [c for c,v in values if v == min([v for c,v in values])]
print(best)

注意：这是原始函数的缩小版本，用于示例目的。

【讨论】：

请注意，这是问题中已经提到的蛮力 try-every-possibility 选项。

问题的症结在于如何在最小化策略下使用scipy.optimize 中的东西返回整数答案。仅仅因为问题中提到了这个概念并不意味着有人会知道使用 brute 或考虑使用 itertools - 尤其是初学者。该步长最初应该为 1 也不是很明显。如果您有更好的答案，请务必发布！

非常感谢 - 肯定会研究纸浆以解决未来的优化问题，不知何故在此之前不知道它存在！

【参考方案2】：

可能对您的问题有所帮助的一件事，您可以有一个约束：

max([x-int(x)])=0

这不会完全解决您的问题，该算法仍会尝试作弊，您将获得具有某种程度错误~±5e-10 的值，它仍会尝试仅通过 scipy 算法中的错误进行优化，但它是总比没有好。

cons = ('type':'eq', 'fun': con,
        'type':'eq','fun': lambda x : max([x[i]-int(x[i]) for i in range(len(x))]))

在我知道解决方案的一些优化上测试了这个过程，这个过程对初始值比无约束搜索更敏感，它得到相当准确的答案但是解决方案实际上可能找不到真实值，你基本上需要优化过程的大跳跃（它用来确保它没有优化到局部最小值）来搜索样本空间，因为较小的增量通常不足以移动到下一个数字。

【讨论】：

有趣的想法——正如你所说，不是一个完整的解决方案，但总比没有好，谢谢！

【参考方案3】：

以下是使用Python Gekko（我维护的一个包）解决混合整数非线性规划问题的方法：

from gekko import GEKKO

m = GEKKO(remote=False)
x = m.Array(m.Var,9,lb=0,ub=7,integer=True)

def f(x):
    return (481.79/(5+x[0]))+(412.04/(4+x[1]))\
           +(365.54/(3+x[2]))+(375.88/(3+x[3]))\
           +(379.75/(3+x[4]))+(632.92/(5+x[5]))\
           +(127.89/(1+x[6]))+(835.71/(6+x[7]))\
           +(200.21/(1+x[8]))

m.Minimize(f(x))
m.Equation(sum(x)==7)
m.options.SOLVER=1
m.solve()
print(x)

这给出了解决方案：

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.0529 sec
 Objective      :  859.5269999999999
 Successful solution
 ---------------------------------------------------


[[0.0] [0.0] [1.0] [1.0] [1.0] [0.0] [1.0] [0.0] [3.0]]

【参考方案4】：

暴力破解的通用函数。在 scipy 中比 brute 做得更好，因为 scipy 实际上使用浮点数运行函数，而不仅仅是整数，尽管 range 明确说明了这一点，正如 Jarad 所述

def brute(func, arg_ranges, finish=min):
if isinstance(arg_ranges, dict):
    args = k:np.unique(np.hstack([a for r in rs for a in r]) if isinstance(rs, list) else [a for a in rs]) for k,rs in arg_ranges.items()
    print(args)
    return finish([(dict(zip(args.keys(), vs)), func(**dict(zip(args.keys(), vs)))) for vs in itertools.product(*args.values())], key=lambda x: x[1])
elif isinstance(arg_ranges, list):
    return finish([(i, func(i)) for r in arg_ranges for i in r], key=lambda x: x[1])
else:
    return finish([(i, func(i)) for i in arg_ranges], key=lambda x: x[1])

print(brute(lambda x,y: x / (y + 2), 'x':[range(1,5,2), range(0,6,1)], 'y':range(2,5,1)))