实施具有自定义约束的 Python 求解器

Question

我有两个相互关联的变量，我想找到一个最优解，在本例中是它们总和的最小值。现在，我们称它们为 X 和 Y，加上预定义的常量，它们加起来就是一组 "variables" s1 和 s2（稍后提供约束）：

105896649.59 + X = s1
    -6738.82 + Y = s2

在搜索 SciPy 文档时，我遇到了 linear programming 解决方案，其中我有最小化函数（在本例中为 X + Y）和一组不等式和我的变量所绑定的等式约束。就我而言，它们如下：

X >= 0, Y >= 0
s1 >= 1, s2 >= 1
s2 / (s1 + s2) = 0.0001%

对于这个具体案例，代码很容易实现：

from scipy.optimize import linprog

lstConst = [105896649.59, -6738.82]

# function to minimise: X + Y
c= [1, 1]

# left-hand side of the equation for s2 / (s1 + s2) = 0.0001%
# i.e., -0.000001 * X + 0.999999 * Y
Aeq = [[-0.000001, 0.999999]]

# right-hand side of the equation
beq = [0.000001 * (lstConst[0] + lstConst[1]) - lstConst[1]]

# ensures minimum can't be a negative number
minX = max(1, max(1 -lstConst[0], 0))
minY = max(1, max(1 -lstConst[1], 0))

X_bounds = (minX, None)

Y_bounds = (minY, None)

res = linprog(c, A_eq=Aeq, b_eq=beq, bounds=[X_bounds, Y_bounds])

所以我们有 X 和 Y 的值来最小化 x 参数上的函数：

In [1]: res.x
Out[1]: array([1.00000000e+00, 6.84471676e+03])

我想以此方法为基础：

其实还有另外一组限制：s1和s2也必须是整数（注意X和Y是浮点数没有问题).
我不会为 s1 和 s2 之间的比率定义单个值，而是提供一个包含不同可能比率的列表。

本质上，我想在给定 s1 和 s2 之间的几个不同比率的情况下找到 X + Y 函数的最小值。这可以通过遍历列表以在每次迭代中定义 Aeq 和 beq 或定义其他限制（如果可能）来实现。

但是，我对整数限制以及如何使线性规划算法考虑到它一无所知。

如果有人有使用 library/optimizer 而不是 SciPy 和 linprog 的替代建议，我们也欢迎。

Answer 1

首先重申问题：

minimize x + y, subject to:

    k1 + x = s1
    k2 + y = s2
    x >= 0
    y >= 0
    s1 >= 1
    s2 >= 1
    s2 / (s1 + s2) = k3

Where:

    k1 = 105896649.59
    k2 = -6738.82
    k3 = 0.000001

请注意，您不需要 s1 和 s2 变量来对 linprog 中的问题进行编码。如果没有 s1 和 s2 辅助变量，问题是：

minimize x + y, subject to:

  x >= 0
  y >= 0
  x + k1 >= 1,
  y + k2 >= 1,
  (1-k3)y - k3x = (k1 + k2)k3 - k2

在 linprog 中更容易阅读和编码：

import numpy as np
from scipy.optimize import linprog
k1, k2, k3 = 105896649.59, -6738.82, 0.000001
A_ub = -np.eye(2)
b_ub = [k1-1, k2-1]
A_eq = [[-k3, (1-k3)]]
b_eq = (k1 + k2)*k3 -k2
res = linprog([1,1], A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=[[0,None], [0, None]])
print(res)

这给出了 [0., 6844.71675549] x=1 的位置，因为您实际上已经将 x 和 y 的下限设置为 1（我认为这是一个错字......）但它没有在所问问题的上下文中无关紧要：

进入问题：

... I'm clueless as to the integer restriction and how to make the linear programming algorithm take it into account.

If anyone has an alternative suggestion that uses a library/optimizer other than SciPy and linprog, that's also welcome.

你要的是mixed integer linear programming (MILP). MILP and linear programming (LP), are typically solved with different algorithms, and a MILP problem is typically harder to solve exactly. SciPy Optimize doesn't support MILP. There are a number of open source tools that do such as OrTools and PySCIPOpt which is a Python wrapper over SCIP.

PySCIPOpt 中的示例：

PySCIPOpt 很好，因为它有一个约束编程类型 API。在 PySCIPOpt 中，您的问题很容易以可读的形式陈述。重新引入辅助变量，我们几乎可以逐字输入约束：

from pyscipopt import Model

k1, k2, k3 = 105896649.59, -6738.82, 0.000001
model = Model()
x = model.addVar(vtype="CONTINUOUS", name="x", lb=0)
y = model.addVar(vtype="CONTINUOUS", name="y", lb=0)
s1 = model.addVar(vtype="CONTINUOUS", name="s1", lb=None, ub=None)
s2 = model.addVar(vtype="CONTINUOUS" name="s2", lb=None, ub=None)
o = model.addVar(vtype="CONTINUOUS", name="Objective Value", lb=0, ub=None)
model.addCons(k1 + x == s1)
model.addCons(k2 + y == s2)
model.addCons(s1 >= 1)
model.addCons(s2 >= 1)
model.addCons(s2/(s1+s2) == k3)
model.addCons(x + y == o)
model.setObjective(o, "minimize")
model.optimize()
print('x + y = o -> (%.4f + %.4f = %.4f)' % (model.getVal(x), model.getVal(y), model.getVal(o)))

给出与 linprog 相同的答案，因为它只是一个线性程序。然而，由于 SCIP 支持 MILP，我们可以引入整数变量。要处理您的案例 #1，只需将 s1 和 s2 更改为整数：

...
s1 = model.addVar(vtype="INTEGER", name="s1", lb=None, ub=None)
s2 = model.addVar(vtype="INTEGER", name="s2", lb=None, ub=None)

给出：

...
SCIP Status        : problem is solved [optimal solution found]
Solving Time (sec) : 0.00
Solving Nodes      : 1
Primal Bound       : +1.10089229999989e+05 (1 solutions)
Dual Bound         : +1.10089229999989e+05
Gap                : 0.00 %
x + y = o -> (103244.4100 + 6844.8200 = 110089.2300)

这是一个完全不同的解决方案...但这就是 MILP 不是 LP 的原因。

从上面的示例中，通过阅读 docs，您应该能够弄清楚如何编写您的 #2 案例 - 基本上像 1/k3 这样的东西成为您模型中的另一个整数变量。

实施具有自定义约束的 Python 求解器

Implementing a Python Solver with Customised Constraints

python

scipy

linear-programming