解决两个变量线性不等式的算法

Question

我正在尝试找到一种算法，以确定带有两个变量和以下形式的线性不等式集的严格正整数解的存在：

• ₁ + ₁ ≤ ₁
• ₂ + ₂ ≤ ₂
• ₃ + ₃ ≤ ₃
• ...

该问题还涉及以下形式的最终不等式：

• + ≥

在这里，一些线性规划技术应该适用，而我对它们不是很熟悉。我一直在寻找一种更为临时的解决方案，或许可以利用问题中给出的性质和约束（不等式类型）。任何见解或算法都将受到欢迎，但是特别有趣的是一些确定性的和线性的算法，其中 是不等式数量。

附加约束条件： ，， > 0

- Futuristic Gladiator

---

如果只有两个变量，那么我认为你只需要在所有地方使用 x 和 y，而不是 x1、x2 等等。并且可以说，涉及的所有数字（ai、bi、li、x、y）都大于 0，这样安全吗？ - user3386109

哦，你说得对，抱歉。虽然是的，所有的ai、bi、li都大于0。 - Futuristic Gladiator

好的，那么x和y可以是任何整数，包括0和负整数？ - user3386109

从技术上讲，但请注意最终的不等式和对li的限制，因为意味着解决方案必须大于0，但是它在技术上可以是任何整数。 - Futuristic Gladiator

可能有无限的解决方案。是否有一些表达式可以进行优化（最大化、最小化）？ - trincot

Answer 1

第一组不等式描述了斜率为负的直线。换句话说，它们的形式为=+，其中是一个负有理数（这来自于，， > 0），它们排除了其上方的区域。

最后一个不等式描述了一条斜率为负45°的直线，它排除了其下方的区域。

我们可以像这样描绘一个例子问题：

白色区域包含解点。黑色区域被第一组不等式排除，紫色区域被最后一个不等式排除，而紫色区域则排除了负或负的解。

由于橙色线的斜率保证为负，因此我们可以得出结论：如果最后一条线（斜率为-1的深红色线）上没有属于解的点，则也没有其他点。没有其他点能够在解中，同时最后一条（深红色）线上的点都被排除。

所以，我们“仅需”检查最后一行上的整数点。换句话说，我们可以用以下约束替换最终约束:

      + =

我们可以遍历所有其他约束，并查看相应线的斜率与最后一行的斜率如何比较：

• 如果小于：橙色线比最后一行更“水平”。与最后一行的交叉表示的是的下限(和 的上限)。
• 如果大于：橙色线比最后一行更“垂直”。与最后一行的交叉表示的是 的上限(和 的下限)。
• 如果相等：橙色线与最后一行平行。如果它在最后一行上方(或与其相等)，则可以忽略它。如果它在最后一行下方，则没有解。

在此过程结束时，我们有的最大下限和最小上限。如果这些(有理数)边界仍然允许正整数 的解，并且相应的 也是正数，则报告成功。否则，就没有解。

该算法的时间复杂度为O()。

- trincot

Answer 2

为了完整起见，我展示一个传统的线性规划解决方案。优点是它相当容易构建。缺点是，虽然它会高效，但它不太可能像在问题中寻找手动、分析或半分析解决方案那样高效。

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


def solve(ab: np.ndarray, l: np.ndarray) -> OptimizeResult:
    lb = np.full(n+1, -np.inf)  # No lower bounds...
    lb[-1] = l[-1]              # except for the last constraint
    ub = np.full_like(lb, np.inf)  # Default: no upper bounds (for last constraint)
    ub[:-1] = l                    # Most upper bounds determined by 'l'
    A = np.ones((n+1, 2))  # For last constraint: coefficients are 1
    A[:-1, :] = ab         # All other coefficients are 'ab'

    return milp(
        c=np.zeros(2),                 # no optimisation, only bounds-testing
        integrality=np.ones(2),        # x, y both integral
        bounds=Bounds(lb=1),           # only a lower simple bound, of 1
        constraints=LinearConstraint(lb=lb, A=A, ub=ub),
    )


def explain(ab: np.ndarray, l: np.ndarray) -> None:
    a, b = ab.T
    print('a =', a)
    print('b =', b)
    print('l =', l)
    result = solve(ab, l)
    print(result.message)
    if result.success:
        print('x, y =', result.x)
        print('[a b][x y] =', ab @ result.x)
        print('          <=', l)
        print(f'x + y = {result.x.sum()} >= {l[-1]}')
    print()


# Infeasible example
rand = default_rng(seed=1)
n = 6
abl = rand.uniform(low=0.1, high=1, size=(n, 3))
ab = abl[:, :2]
explain(ab, l=abl[:, 2])

# Feasible example
explain(ab, l=ab @ (10, 10))

a = [0.56063946 0.9537845  0.84493233 0.1248032  0.39675854 0.5081481 ]
b = [0.95541733 0.38064831 0.46827922 0.7781618  0.80958583 0.22063753]
l = [0.22974365 0.4809938  0.59463432 0.58432898 0.37287535 0.46280169]
The problem is infeasible. (HiGHS Status 8: model_status is Infeasible; primal_status is At lower/fixed bound)

a = [0.56063946 0.9537845  0.84493233 0.1248032  0.39675854 0.5081481 ]
b = [0.95541733 0.38064831 0.46827922 0.7781618  0.80958583 0.22063753]
l = [15.16056789 13.34432809 13.13211557  9.02965    12.06344378  7.28785628]
Optimization terminated successfully. (HiGHS Status 7: Optimal)
x, y = [7. 1.]
[a b][x y] = [4.87989356 7.05713982 6.38280556 1.65178421 3.58689565 3.77767423]
          <= [15.16056789 13.34432809 13.13211557  9.02965    12.06344378  7.28785628]
x + y = 8.0 >= 7.287856280550347

- Reinderien