Solving an integer (boolean) constraint satisfaction problem I have a 0-1 integer constraint satisfaction problem of the following form: find binary vectors $x = (x_1,\dots,x_m) \in \{0,1\}^m$ and $y = (y_1, \dots,y_n) \in \{0,1\}^n$ that satisfy the constraints


*

*$x_i \le \sum_{j,k} a_{ijk} x_j y_k\ $ for $i = 1,\dots,m$

*$x_i \ge a_{ijk} x_j y_k\ $ for $i = 1,\dots,m$, $j = 1,\dots,m$, $k = 1,\dots,n$

*$\sum_j b_{lj} x_i \ge c_l\ $ for $l = 1,\dots,p$
where $a_{ijk}$, $b_{lj}$ and $c_l$ are constants, known beforehand. The $a_{ijk}$ also take values in $\{0,1\}$, and most of their elements will be zero (they are sparse arrays). However, $b_{lj}$ and $c_l$ can be any positive integers.
Ideally I would want to find all vectors $x$ and $y$ that satisfy the constraints. Is this a known / well understood type of problem? Are there any known methods for solving it, ideally with a solver available?

Note: constraints 1 and 2 are derived from the corresponding boolean constraint
$$x_i = \bigvee_{j,k} a_{ijk} \land x_j \land y_k $$
interpreting the $0/1$ variables as booleans.

The problem you describe is a non-convex binary program. The non-convexity comes from the first and second set of constraints: In $x_i \le \sum_{j,k} a_{ijk} x_j y_k\ $ two decision variables $x_j$ and $y_k$ are multiplied by each other. This will make the problem very hard to solve for most solvers. However, there is a way to turn your problem into a convex problem. 
The product of two binary variables $x$ and $y$ can be linearized by introducing a new variable $z$ with the additional constraints
$$z \le x$$
$$z \le y$$
$$z \ge x+y-1$$
$$0 \le z \le 1$$
Now, if $x$ or $y$ is zero, so is $z$. If both $x=y=1$, then $z=1$. Hence, $z=xy$.
We will apply this approach to your problem now. The first set of constraints contains the product $x_j y_k$. We will introduce a new variable $z_{jk}$ in the above fashion and replace your first set of constraints with
$$x_i \le \sum_{j,k} a_{ijk} z_{jk}$$
$$z_{jk} \le x_j$$
$$z_{jk} \le y_k$$
$$z_{jk} \ge x_j+y_k-1$$
$$0 \le z_{jk} \le 1$$
Doing the same for the second set of constraints will result in the overall problem formulation:
$$x_i \le \sum_{j,k} a_{ijk} z_{jk}$$
$$z_{jk} \le x_j$$
$$z_{jk} \le y_k$$
$$z_{jk} \ge x_j+y_k-1$$
$$0 \le z_{jk} \le 1$$
$$x_i \ge a_{ijk} z_{jk}\ $$
$$ \sum_j b_{lj} x_i \ge c_l\ $$
This formulation can now be plugged into any integer-programming solver. The introduction of the variables $z_{jk}$ has now increased the number of variables in the optimization problem but it turned the initial problem into a convex problem which makes it much easier to solve despite the higher amount of variables.
