用mystic解一个不等式系统

Question

我正在尝试使用Python框架Mystic来解决一个简单的不等式系统。我正在求解十个变量，这些变量都应该大于其他变量(x1 > x2 > x3 > ... > x10)，并且它们应该大于零。我希望这些变量仅限于整数，这是我还没有做到的。

这十个变量被用在一个具有大约15-20个方程的不等式系统中。我已经用所有这些不等式做了一个字符串，简化了它们，生成了求解器，生成了约束，就像文档中解释的那样。然而，我无论如何也不能理解如何从我的代码生成的约束对象中获得解决方案(即，mystic分配给我想要求解的变量的数字)。我已经搜索了几个小时的文档，并检查了大多数关于神秘和不平等解决的stackoverflow帖子。大多数问题都比我的级别高得多，所以它们没有多大帮助。

下面是我的代码摘录。

inequalities = '''
x1 > x2
x2 > x3
# etc.
x9 > x10
x10 > 0
4*x1 + 4*x2 + 2*x3 + 2*x4 + 2*x6 + x7 + 2*x8 > x4 + x6 + x8 + x9
4*x1 + 4*x2 + 2*x3 + 2*x4 + 2*x6 + x7 + 2*x8 > 2*x1 + x2 + 2*x3 + x5 + 2*x6 + x7 + x8 + x10
4*x1 + 4*x2 + 2*x3 + 2*x4 + 2*x6 + x7 + 2*x8 > 2*x1 + 5*x3 + 3*x4 + x8 + x9 + x10
# etc.
'''

def solve_inequalities(inequalities):
    var = ms.get_variables(inequalities)
    eqns = ms.simplify(inequalities, variables=var)
    print("Simplified equations: ", eqns)

    solver = ms.generate_solvers(eqns, var)
    constraint = ms.generate_constraint(solver)
    # which of these objects do i need to call?
    solution = constraint([1,2,3,4,5,6,7,8,9,10])
    print("Solution: ", constraint)

solve_inequalities(inequalities)

打印的内容：

Simplified equations:  x6 > x7
x3 > x4
x2 > -x1 - x3/2 - x4/2 - x6/4 + 3*x9/4 + x10/4
x3 > -2*x1 - 3*x2/2 - x4 + x5 - x6/2 - x7/2 - x8/2 + x9 + 3*x10
x1 > -x2 - x3/4 - x4/2 - x6/2 - x8/4
x2 > -x1 - x3/4 - x4/4 - x6/4 - x8/2 + 3*x9/4 + x10/2
x5 < -2*x1 + x2 - x3 + x4 + x6 + x8 - x10
x5 > x6
x4 > -x1 - 3*x2/2 + x5/2 - x8/2 + x10/2
x4 > x5
x2 > x3
x8 > x9
x5 < 3*x1/2 + x2 + x3/2 + x4/2 - x6/2
x2 > -x1 + x4/4 + x5 - x6/4 - x7/4 - x8/4 + x9/4 + x10/4
x1 > -x2 - x3/2 - x4/4 - x6/2 - x7/4 - x8/2 + x9/2
x1 > -x2 - x3/4 - x4/4 + x5/4 + 3*x7/4
x1 > x2
x2 > -3*x1/4 - x3/2 - x4/2 - x6/2 - x7/4 - x8/2 + x9/2
x1 > -x2 - x3/2 - x4/2 - x6/2 - x7/4 - x8/4
x10 > 0
x7 > x8
x3 > -3*x1/2 - 3*x2/2 - x4/2 + 3*x5/2 - x6 + 3*x7/2 - x8/2 + x10
x3 > -3*x1/2 + x2 - x4/2 + 3*x5/2 - x6/2 + x7/2 - x8/2 + x10
x9 > x10
x5 < 3*x1/2 + 3*x2/2 + x3/2 + x8 - 3*x9/2 - x10/2
x2 > -x1 - x3/2 - x4/4 - x6/4 - x8/4 + x9/4 + x10/4
x2 > -x1/2 + 3*x3/4 + x4/4 - x6/2 - x7/4 - x8/4 + x9/4 + x10/4
x1 > -x2 - x3/2 - x4/4 - x6/4 - x7/4 - x8/4 + x9/4

Solution:  [20.875000000000057, 20.875000000000036, 20.875000000000014, 6.000000000000021, 6.000000000000014, 6.000000000000007, 8.000000000000018, 8.000000000000009, 9.00000000000001, 9, 10]

这显然是不正确的。(我还想知道:为什么这些不等式的原始顺序没有得到保留？)

任何帮助都将不胜感激。

向您致敬，Lennart

Answer 1

我是mystic的作者。很抱歉，文档并不像预期的那样明显。我想，文档中的几个地方有两个微妙的地方，包括generate_constraint。这意味着，如果您有x1 = x2 + x3，然后是x1 = x2 + x4，它将首先将x1设置为x2 + x3，然后将x1设置为x2 + x4……从而替换x1的初始值。要获得您期望的结果，您需要在简化时更改顺序以获得x1 = x2 + x3和x4 = x1 - x2……否则，您需要将join设置为and_。前者可以使用simplify的一些关键字自动尝试，如cycle和target...但不能保证它会比你用sort手动排序公式做得更好。下面我将使用您的简化不等式来演示后者(因为我没有您的原始方程)。

inequalities = '''
x1 > x2
x2 > x3
x3 > x4
x4 > x5
x5 > x6
x6 > x7
x7 > x8
x8 > x9
x9 > x10
x10 > 0
x2 > -x1 - x3/2 - x4/2 - x6/4 + 3*x9/4 + x10/4
x3 > -2*x1 - 3*x2/2 - x4 + x5 - x6/2 - x7/2 - x8/2 + x9 + 3*x10
x1 > -x2 - x3/4 - x4/2 - x6/2 - x8/4
x2 > -x1 - x3/4 - x4/4 - x6/4 - x8/2 + 3*x9/4 + x10/2
x5 < -2*x1 + x2 - x3 + x4 + x6 + x8 - x10
x4 > -x1 - 3*x2/2 + x5/2 - x8/2 + x10/2
x5 < 3*x1/2 + x2 + x3/2 + x4/2 - x6/2
x2 > -x1 + x4/4 + x5 - x6/4 - x7/4 - x8/4 + x9/4 + x10/4
x1 > -x2 - x3/2 - x4/4 - x6/2 - x7/4 - x8/2 + x9/2
x1 > -x2 - x3/4 - x4/4 + x5/4 + 3*x7/4
x2 > -3*x1/4 - x3/2 - x4/2 - x6/2 - x7/4 - x8/2 + x9/2
x1 > -x2 - x3/2 - x4/2 - x6/2 - x7/4 - x8/4
x3 > -3*x1/2 - 3*x2/2 - x4/2 + 3*x5/2 - x6 + 3*x7/2 - x8/2 + x10
x3 > -3*x1/2 + x2 - x4/2 + 3*x5/2 - x6/2 + x7/2 - x8/2 + x10
x5 < 3*x1/2 + 3*x2/2 + x3/2 + x8 - 3*x9/2 - x10/2
x2 > -x1 - x3/2 - x4/4 - x6/4 - x8/4 + x9/4 + x10/4
x2 > -x1/2 + 3*x3/4 + x4/4 - x6/2 - x7/4 - x8/4 + x9/4 + x10/4
x1 > -x2 - x3/2 - x4/4 - x6/4 - x7/4 - x8/4 + x9/4
'''

因此，就像您之前所做的一样，但这一次，使用join=and_。

>>> import mystic as my
>>> var = my.symbolic.get_variables(inequalities)
>>> constraint = my.symbolic.generate_constraint(my.symbolic.generate_solvers(inequalities, var), join=my.constraints.and_)
>>> constraint([1,2,3,4,5,6,7,8,9,10])
[12.750000000000028, 12.750000000000014, 12.75, 10.000000000000064, 10.000000000000053, 10.000000000000043, 10.000000000000032, 10.000000000000021, 10.00000000000001, 10]

你可以看到它起作用了。

如果你想构建一个尝试使用整数的约束，那么你可以这样做：

>>> c = my.constraints.integers(float)(lambda x:x)
>>> c([1.1,2.3,3.7,4.3,5,6,7,8.932,9.0002,10])
[1.0, 2.0, 4.0, 4.0, 5.0, 6.0, 7.0, 9.0, 9.0, 10.0]

然而，将这两个约束耦合在一起可能会失败，正如您在这里看到的：

>>> con = my.constraints.and_(c, constraint, maxiter=1000)
>>> con([1,2,3,4,5,6,7,8,9,10])
[12.750000000000028, 12.750000000000014, 12.75, 10.000000000000064, 10.000000000000053, 10.000000000000043, 10.000000000000032, 10.000000000000021, 10.00000000000001, 10.0]

它失败的原因是，这是另一个微妙之处，因为and_本质上试图循环使用迭代应用约束的组合……并希望它能成功。通常是这样的，但在你的情况下，它不会...因此，它放弃并返回原始的x。

事实上，即使这样也不是成功的：

>>> ineq = '''
... x1 > x2
... x2 > x3
... x3 > x4
... x4 > x5
... x5 > x6
... x6 > x7
... x7 > x8
... x8 > x9
... x9 > x10
... x10 > 0
... '''
>>> 
>>> constraint = my.symbolic.generate_constraint(my.symbolic.generate_solvers(ineq, var), join=my.constraints.and_)
>>> con = my.constraints.and_(c, constraint, maxiter=1000)
>>> con([1,2,3,4,5,6,7,8,9,10])
[10.000000000000096, 10.000000000000085, 10.000000000000075, 10.000000000000064, 10.000000000000053, 10.000000000000043, 10.000000000000032, 10.000000000000021, 10.00000000000001, 10.0]
>>> 

因此，知道了这一点，我们可以像这样重写约束：

>>> ineq = '''
... x1 > 0
... x1 < x2
... x2 < x3
... x3 < x4
... x4 < x5
... x5 < x6
... x6 < x7
... x7 < x8
... x8 < x9
... x9 < x10
... '''
>>> constraint = my.symbolic.generate_constraint(my.symbolic.generate_solvers(ineq, var), join=my.constraints.and_)
>>> con = my.constraints.and_(c, constraint, maxiter=10)
>>> con([1,2,3,4,5,6,7,8,9,10])
[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
>>> con([1.1,2.3,3.7,4.3,5,6,7,8.932,9.0002,10])
[1.0, 2.0, 3.999999999999995, 4.0, 5.0, 6.0, 7.0, 8.99999999999999, 9.0, 10.0]

你可以看到它做得更好，而且这一次几乎可以工作了。但在第二种情况下仍然失败。

所以，在这种情况下...如果你打算在优化问题中使用这些约束，我会使用constraint作为硬约束(如上所述)，然后应用整数约束作为惩罚(例如软约束)……反之亦然。当我有两个and_约束失败时，我通常会应用其中一个作为惩罚。

如果mystic能更好地应用来自不平等系统的约束，再加上额外的约束(比如整数约束)，那就更好了，但它主要是在幕后蛮力，所以它可能会失败，迫使您依赖于使用penalty。

如果您有合理的界限，对于x，您可以执行以下操作...

>>> cu = my.constraints.and_(my.constraints.discrete(range(100))(lambda x:x), my.constraints.impose_unique(range(100))(lambda x:x))
>>> 
>>> cu([1.1,2.3,3.9,4.3,5,6,7,8.9,9.0002,10])
[1.0, 2.0, 4.0, 38.0, 5.0, 6.0, 7.0, 9.0, 75.0, 10.0]
>>> 
>>> con = my.constraints.and_(c, constraint, cu, maxiter=1000)
>>> con([1.1,2.3,3.9,4.3,5,6,7,8.9,9.0002,10])
[1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 9.0, 10.0, 32.0, 47.0]
>>> 

因此，这意味着，从0,99中的离散整数集中选择解决方案，其中所有项目都是唯一的……

所以，现在，回到你的更难的原始案例...

>>> constraint = my.symbolic.generate_constraint(my.symbolic.generate_solvers(inequalities, var), join=my.constraints.and_)
>>> con = my.constraints.and_(c, constraint, cu, maxiter=1000)
>>> con([1.1,2.3,3.9,4.3,5,6,7,8.9,9.0002,10])
[99.0, 98.0, 97.0, 96.0, 67.0, 54.0, 29.0, 20.0, 6.0, 1.0]

>>> x1, x2, x3, x4, x5, x6, x7, x8, x9, x10 = [99.0, 98.0, 97.0, 96.0, 67.0, 54.0, 29.0, 20.0, 6.0, 1.0]
>>> x1 > -x2 - x3/2 - x4/4 - x6/4 - x7/4 - x8/4 + x9/4
True
>>> x2 > -x1/2 + 3*x3/4 + x4/4 - x6/2 - x7/4 - x8/4 + x9/4 + x10/4
True
>>> x2 > -x1 - x3/2 - x4/4 - x6/4 - x8/4 + x9/4 + x10/4
True

耶！

编辑:有趣的是，通过构建一个小函数来检查哪个不等式失败了…

def failures(x):
  return tuple(eval(i, dict(zip(var,x))) for i in inequalities.strip().split('\n'))

x5 < -2*x1 + x2 - x3 + x4 + x6 + x8 - x10

当包含此行时，不等式求解器运行在接近最大1000次尝试的情况下...然而，如果我取消了这个限制，那么成功的概率是100%，而且非常快(求解器使用了很少的尝试)。

因此，这里的结论是应该删除这一行，并将其重铸为惩罚…而其他一切都应该作为约束来处理。