给定一个线性系统Ax = b
,其中矩阵A
和向量b
具有整数值,我想找到所有非负整数向量x
求解这个方程。
到目前为止,我已经发现了一些技术,例如史密斯范式或矩阵的埃尔米特范式来找到整数解,我想我可以使用线性求解器来找到非负解。有没有一个库可以让这更容易?
Python 解决方案将是理想的,但如果一个库存在于另一种语言中,我想知道它。
您可以使用整数规划来执行此操作,为每个 x 值定义一个非负整数决策变量。然后,您可以使用约束 Ax=b 和目标 0 来解决问题,这会搜索方程组的任何可行整数解。
这可以使用 python 中的 pulp 包轻松实现。例如,考虑求解以下系统:
x+2y+z = 5
3x+y-z = 5
您可以通过以下方式在纸浆中解决此问题:
import pulp
A = [[1, 2, 1], [3, 1, -1]]
b = [5, 5]
mod = pulp.LpProblem('prob')
vars = pulp.LpVariable.dicts('x', range(len(A[0])), lowBound=0, cat='Integer')
for row, rhs in zip(A, b):
mod += sum([row[i]*vars[i] for i in range(len(row))]) == rhs
mod.solve()
[vars[i].value() for i in range(len(A[0]))]
# [1.0, 2.0, 0.0]
这标识了一个可行的整数解,x=1,y=2,z=0。
如果至少存在一个Ax=b
的非负积分解,那么整数规划应该能够找到它。所有非负积分解的集合可以通过A
的零空间找到。
示例
使用欧文回答中的A
和b
:
>>> from sympy import *
>>> A = Matrix([[ 1, 2, 1],
[ 3, 1,-1]])
>>> b = Matrix([20,12])
计算空空间:
>>> A.nullspace()
[Matrix([
[ 3/5],
[-4/5],
[ 1]])]
零空间是一维的,其基向量是有理值的。使用 Erwin 找到的特定解(2,8,2)
,所有实解的集合是一条参数化为(2,8,2) + t (3,-4,5)
的线。由于我们只对非负积分解感兴趣,因此我们将直线与非负辛烷相交,然后与整数格相交,从而产生t ∈ {0,1,2}
。因此,如果:
t=1
,我们得到(5,4,7)
的解。t=2
,我们得到(8,0,12)
的解。
请注意,这些是 Erwin 使用 Z3 找到的 3 种解决方案。
但是,当零空间的维数大于 1 时,这要困难得多。看看:
- Jesús De Loera,多面体中计算晶格点的许多方面,2005年。
这是一个 Z3 模型。Z3 是来自 Microsoft 的定理证明器。该模型与之前提出的MIP模型非常相似。
在 MIP 中枚举整数解并非完全简单(可以通过一些努力 [link] 或使用高级 MIP 求解器中的"解决方案池"技术来完成)。使用 Z3,这要容易一些。使用约束规划(CP)求解器可能更容易:它们可以自动枚举解决方案。
来吧:
from z3 import *
A = [[1, 2, 1], [3, 1, -1]]
b = [20, 12]
n = len(A[0]) # number of variables
m = len(b) # number of constraints
X = [ Int('x%d' % i) for i in range(n) ]
s = Solver()
s.add(And([ X[i] >= 0 for i in range(n) ]))
for i in range(m):
s.add( Sum([ A[i][j]*X[j] for j in range(n) ]) == b[i])
while s.check() == sat:
print(s.model())
sol = [s.model().evaluate(X[i]) for i in range(n)]
forbid = Or([X[i] != sol[i] for i in range(n)])
s.add(forbid)
它解决了
x0+2x1+x2 = 20
3x0+x1-x2 = 12
解决方案如下所示:
[x2 = 2, x0 = 2, x1 = 8]
[x2 = 7, x1 = 4, x0 = 5]
[x2 = 12, x1 = 0, x0 = 8]
我们可以打印最终模型,以便我们可以看到"禁止"约束是如何添加的:
[And(x0 >= 0, x1 >= 0, x2 >= 0),
1*x0 + 2*x1 + 1*x2 == 20,
3*x0 + 1*x1 + -1*x2 == 12,
Or(x0 != 2, x1 != 8, x2 != 2),
Or(x0 != 5, x1 != 4, x2 != 7),
Or(x0 != 8, x1 != 0, x2 != 12)]
请参阅此处zsolve
4ti2
包的方法。