Using mystic to solve a system of inequalities

Question

I am trying to solve a system of simple inequalities using the Python framework Mystic. I'm solving for ten variables, these variables should all be greater than each other (x1 > x2 > x3 > ... > x10) and they should be greater than zero. I'd like for these variables to be restricted to only integers, something which I've not managed to get done yet.

These ten variables are used in a system of inequalities with about 15-20 equations. I've made a string with all these inequalities, simplified them, generated solvers, generated constraints out of these solvers, just like the documentation explains. However, I cannot for the life of me understand how to get the solution (that is, the numbers mystic assigns to the variables for which I want to solve) out of this constraint object my code generates. I've scoured the documentation for hours now, and checked most of the stackoverflow posts regarding mystic and inequality solving. Most questions are of a much higher level than mine, so they aren't of much help.

Here's an excerpt of my code.

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)

Which prints:

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]

Which is obviously not correct. (Something else I'm wondering: Why does the original order of these inequalities not get preserved?)

Any help would be greatly appreciated.

Regards, Lennart

Answer 1

I'm the mystic author. Sorry that the docs are not obvious as expected. There's, I guess, two subtle points that are in the documentation in several places, including in generate_constraint. One is that the constraints are applied iteratively, by default, and it's the responsibility of the user to ensure that the individual equations don't conflict. What that means, is that if you have x1 = x2 + x3 and then x1 = x2 + x4, it will first set x1 to x2 + x3 then set x1 to x2 + x4... thus replacing the initial value of x1. To get what you expected, you'd need to either change the order when you simplify to obtain x1 = x2 + x3 and x4 = x1 - x2... or you'd need to set join to and_. The former can be attempted automatically with some of the keywords for simplify, like cycle and target... but there's no guarantee it will do better than you manually ordering the equations with a sort. I'll demonstrate the latter below, using your simplified inequalities (since I don't have your original equations).

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
'''

So, just as you did before, but this time, using 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]

You can see it works.

If you'd like to build a constraint that attempts to use integers, then you can do this:

>>> 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]

However, coupling these two constraints together is going to likely fail, as you can see here:

>>> 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]

The reason it fails, which is the other subtle point, is that and_ essentially tries to cycle through combinations of iteratively applying the constraints... and hopes that it succeeds. Often it does, but in your case, it doesn't... so it gives up and returns the original x.

Matter of fact, not even this is successful:

>>> 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]
>>> 

However, you'll note that it seems that mystic is keying on the 10 at the end of the constraints and then making sure that each entry before it is slightly larger. So, knowing that, we can rewrite the constraints like this:

>>> 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]

You can see that it does much better, and nearly works this time... but still fails in the second case.

So, in this case... if you are going to be using these constraints in an optimization problem, I'd use constraint as a hard constraint (as built above), and then apply the integer constraint as a penalty (e.g. a soft constraint)... or vice-versa. When I have two constraints that fail with and_, I typically apply one as a penalty.

It would be nice if mystic was better at applying constraints from systems of inequities coupled to additional constraints (like integer constraints), but it's mostly brute force under the covers, so it can fail and force you to rely on using a penalty.

EDIT: Aha... I had a second thought. If you have reasonable bounds, for the x, you can do the following... which has a better shot at working.

>>> 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]
>>> 

So, this says, pick solutions from the discrete set of integers in [0,99], where all the items are unique... and make sure they satisfy your inequality constraints.

So, now, going back to your much harder original case... it might work:

>>> 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]

and checking... it looks good so far...

>>> 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

Yay!

EDIT: Interestingly, by building a little function to check which of the inequalities fails... I noted something interesting.

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

Essentially yields that the solutions sometimes failed, even with the new constraints solver -- but the failure was only in one equation, and only sometimes.

This one is apparently difficult:

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

When this line is included, the inequality solver runs near the maximum 1000 tries... sometimes succeeding, and sometimes failing. However, if I remove this constraint, then success happens 100% of the time, and very quickly (with the solver using few tries).

So, the conclusion here is that this line should be removed, and recast as a penalty... and everything else should be handled as a constraint.