是否可以使用Python混合整数线性规划来获取区间内的所有解决方案？

Question

我有一个线性问题来解决寻找整数。我找到了一种使用 spicy 中新的 milp 实现来解决它的方法。以下是演示代码。

问题如下。从权重向量 w 我正在寻找整数向量 x 例如 x 和权重的点积在给定范围内。看起来像这样

# minimize
abs(w^T @ x - target)

我将其翻译成以下内容以在 milp 中实现：

# maximize
w^T @ x
# constraints
target - error <= w^T @ x <= target + error

在我的特定上下文中，x 可能存在多种解决方案。有没有办法在给定的时间间隔内获得所有解决方案而不是最大化（或最小化）某些东西？

这是 milp 的实现。

import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds

# inputs
ratio_l, ratio_u = 0.2, 3.0
max_bounds = [100, 200, 2, 20, 2]
target = 380.2772 # 338.34175
lambda_parameter = 2
error = lambda_parameter * 1e-6 * target

# coefficients of the linear objective function
w = np.array([12.0, 1.007825, 14.003074, 15.994915, 22.989769], dtype=np.float64)

# the aim is to minimize
#    w^T x - target_mass

# instead I maximize
#    w^T x
# in the constraint domain
#    target - error <= w^T x <= target + error

# constraints on variables 0 and 1:
# ratio_l <= x[1] / x[0] <= ratio_u
# translation =>
#    (ratio_l - ratio_u) * x[1] <= -ratio_u * x[0] + x[1] <= 0
#    use max (x[1]) to have a constant

# linear objective function
c = w

# integrality of the decision variables
# 3 is semi-integer = within bounds or 0
integrality = 3 * np.ones_like(w)  

# Matrice A that define the constraints
A = np.array([
    # boundaries of the mass defined from lambda_parameters 
    w,
    # c[1] / c[0]  max value
    [-ratio_u, 1.0, 0., 0., 0.],
])

# b_up and b_low vectors
# b_low <= A @ x <= b_up
n_max_C = max_bounds[0]
b_up = [
    target + error,  # mass target
    0.,   # c[1] / c[0] constraints up
]
b_low = [
    target - error,  # mass target
    (ratio_l - ratio_u) * max_bounds[0],  # H_C constraints up
]

# set up linear constraints
constraints = LinearConstraint(A, b_low, b_up)

bounds = Bounds(
    lb=[0, 0, 0, 0, 0],
    ub=max_bounds,
)

results = milp(
    c=c,
    constraints=constraints,
    integrality=integrality,
    bounds=bounds,
    options=dict(),
)

print(results)

结果是这样的

            fun: 380.277405
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 380.27643944560145
        mip_gap: 2.5390790665913637e-06
 mip_node_count: 55
         status: 0
        success: True
              x: array([19., 40.,  0.,  7.,  0.])

但它存在其他可能的 x 数组但错误率最高。这个是

m = np.dot(w, [19., 40.,  0.,  7.,  0.])
print(f"{'target':>10s} {'calc m':>27s} {'deviation':>27s} {'error':>12s}      match?")
print(f"{target:10.6f} {target - error:14.6f} <= {m:10.6f} <= {target + error:10.6f}"
      f" {m - target:12.6f} {error:12.6f}   -> {target - error <= m <= target + error}")

    target                      calc m                   deviation        error      match?
380.277200     380.276439 <= 380.277405 <= 380.277961     0.000205     0.000761   -> True

这两个其他示例也有效，我想知道如何在不实现网格算法（如 scipy 中的 brute）的情况下获得它们。

m = np.dot(w, [20., 39.,  1.,  4.,  1.])
print(f"{'target':>10s} {'calc m':>27s} {'deviation':>27s} {'error':>12s}      match?")
print(f"{target:10.6f} {target - error:14.6f} <= {m:10.6f} <= {target + error:10.6f}"
      f" {m - target:12.6f} {error:12.6f}   -> {target - error <= m <= target + error}")

    target                      calc m                   deviation        error      match?
380.277200     380.276439 <= 380.277678 <= 380.277961     0.000478     0.000761   -> True

m = np.dot(w, [21., 38.,  2.,  1.,  2.])
print(f"{'target':>10s} {'calc m':>27s} {'deviation':>27s} {'error':>12s}      match?")
print(f"{target:10.6f} {target - error:14.6f} <= {m:10.6f} <= {target + error:10.6f}"
      f" {m - target:12.6f} {error:12.6f}   -> {target - error <= m <= target + error}")

    target                      calc m                   deviation        error      match?
380.277200     380.276439 <= 380.277951 <= 380.277961     0.000751     0.000761   -> True

Answer 1

线性规划用于优化，通常是在解空间中选择单个最佳解。您的问题无法选择单一的解决方案。由于您具有边界明确的整数变量，因此“蛮力”（尽管不是迭代蛮力）非常实用。这看起来像：

• 在 x 的已知范围内，尝试所有维度的所有值，不包括最大的维度（x1，范围在 0 到 200 之间）
• 根据您的“比率”约束计算边界数组
• 根据“来自目标的错误”约束计算边界数组
• 将两者结合起来找到整体边界
• 过滤边界内的积分解

import numpy as np

ratio_l, ratio_u = 0.2, 3.0
max_bounds = (100, 200, 2, 20, 2)
target = 380.2772
lambda_parameter = 2
error = 1e-6 * lambda_parameter * target
w = np.array((12.0, 1.007825, 14.003074, 15.994915, 22.989769))
i0234 = [0, 2, 3, 4]
w0234 = w[i0234]

x0234 = np.stack(np.meshgrid(
    *(np.arange(1+max_bounds[m]) for m in i0234)
)).reshape((4, -1))
x0, x2, x3, x4 = x0234

x1_ratio_lower, x1_ratio_upper = np.multiply.outer((ratio_l, ratio_u), x0)
x1_target_lower, x1_target_upper = (target - np.add.outer((error, -error), w0234@x0234))/w[1]

x1_lower = np.ceil(np.max((x1_ratio_lower, x1_target_lower), axis=0)).astype(int)
x1_upper = np.floor(np.min((x1_ratio_upper, x1_target_upper), axis=0)).astype(int)
ok, = (x1_upper >= x1_lower).nonzero()

for i in ok:
    xi = x0234[:, i]
    for x1 in range(x1_lower[i], x1_upper[i]+1):
        x = [xi[0], x1, *xi[1:]]
        target_approx = w.dot(x)
        error_approx = target_approx - target
        print(f'x={x} w@x={target_approx:.6f} ~ {target}, '
              f'error={error_approx:.2e}<{error:.2e}')

x=[19, 40, 0, 7, 0] w@x=380.277405 ~ 380.2772, error=2.05e-04<7.61e-04
x=[20, 39, 1, 4, 1] w@x=380.277678 ~ 380.2772, error=4.78e-04<7.61e-04
x=[21, 38, 2, 1, 2] w@x=380.277951 ~ 380.2772, error=7.51e-04<7.61e-04