Finding solutions to linear system of equations with integer constraint in scipy

Question

I have a system of equations where each equation is a linear equation with boolean constraints. For example:

x1 + x2 + x3 = 2
x1 + x4 = 1
x2 + x1 = 1

And each x_i is either 0 or 1. Sometimes there might be a small positive (<5) coefficient (for example x1 + 2 * x3 + x4 = 3. Basically a standard linear programming task. What I need to do is to find all x_i which are guaranteed to be 0 and all x_j which are guaranteed to be 1. Sorry if my terminology is not correct here but by guaranteed I mean that if you generate all possible solutions you in all of them all x_i will be 0 and in all of them x_j will be 1.

For example my equation has only 2 solutions:

• 1, 0, 1, 0
• 0, 1, 1, 1

So you do not have guaranteed 0 and have x_3 as a guaranteed 1.

I know how to solve this problem with or-tools by generating all solutions and it works for my usecases (equations are pretty constrained so usually there are < 500 solutions although the number of variables is big enough to make the whole combinatorial search impossible).

The big problem is that I can't use that library (system restrictions above my control) and only libraries available in my case are numpy and scipy. I found that scipy has scipy.optimize.linprog.

It seems like I have found a way to generate one solution

import numpy as np
from scipy.optimize import linprog

A_eq = np.array([
    [1, 1, 1, 0],  # x1 + x2 + x3 = 2
    [1, 0, 0, 1],  # x1 + x4 = 1
    [1, 1, 0, 0]   # x1 + x2 = 1
])
b_eq = np.array([2, 1, 1])
c = np.zeros(4)
bounds = [(0, 1)] * 4

res = linprog(c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs-ipm')
if res.success:
    print(res.x)

But I can't find a way to generate all solutions. Also I am not sure whether there is a better way to do it as all I need to know is to find guaranteed values.

Answer 1

I'm not sure I grasp all the details of your project, but if your goal is to determine which variables are necessarily 0 or 1 (or constants), wouldn't symbolic mathematics help here rather than trying to find many numeric solutions?

For instance, using sympy:

import sympy

x1, x2, x3, x4, x5 = sympy.symbols('x1 x2 x3 x4 x5')
eq1 = sympy.Eq(x1 + x2 + x3, 2)
eq2 = sympy.Eq(x1 + x4, 1)
eq3 = sympy.Eq(x2 + x1, 1)
eq4 = sympy.Eq(x3 + x5, 1)

result = sympy.solve([eq1, eq2, eq3, eq4], (x1, x2, x3, x4, x5))

result is:

{x1: 1 - x4, x2: x4, x3: 1, x5: 0}

Which makes it easy to determine the values that are necessarily 0/1 (or any other constant):

if result:
    for var, expr in result.items():
        if expr.is_constant(): # or: if expr in {0, 1}:
            print(f'{var} is guaranteed to be {expr}')
else:
    print('No possible solution')

Output:

x3 is guaranteed to be 1
x5 is guaranteed to be 0

Alternatively, if you want to input your equations from matrices:

import sympy

A = sympy.Matrix([[1, 1, 1, 0],  # x1 + x2 + x3 = 2
                  [1, 0, 0, 1],  # x1 + x4 = 1
                  [1, 1, 0, 0]   # x1 + x2 = 1
                 ])
b = sympy.Matrix([2, 1, 1])

result = sympy.linsolve((A, b))
# {(1−𝜏0, 𝜏0, 1, 𝜏0)}

constants = {f'x{i}': x for r in result for i, x in
             enumerate(r, start=1) if x.is_constant()}

Output {'x3': 1}