I'm a solving a model that has the following constraint:
$$ c_{p,n} = \sum_{s\in S}\sum_{i \in \{1,2,3\} } x_{p,s,i-1} x_{n,s,i}, \forall (p,n) \in C $$ where both the $c$ and $x$ variables are binary, $S$ is a set of indices. The indices $p$ and $n$ stand for "previous" and "next" and represent the same entity, i.e. they both belong to a set $T$, $p \in T$, $n \in T$. Their combinations are listed in $C$. The constraint above is used to know whether any combination of $p$ and $n$ is chosen. The variables $c_{p,n}$ appear in the cost function weighted by a cost $w_{p,n}\in \mathbb{R}$.
Other relevant constraints in the problem are: $$\begin{align} \sum_{s\in S}\sum_{i \in I} x_{p,s,i} &\leq 1, \forall p \in T \\ \sum_{p\in T} x_{p,s,i} &\leq 1, \forall s \in S, \forall i \in I \\ \sum_{p\in T} x_{p,s,i} &\leq \sum_{p\in T} x_{p,s,i-1}, \forall s \in S, \forall i \in \{1,2,3\}, \\ \end{align}$$ where $I=\{0,1,2,3\}$. Note that in the last constraint, both the left-hand-size and the right-hand-size can sum up to 1 at most (due to the second constraint).
The bilinear constraints above make my problem slow to solve in certain (large) instances. I tried to reformulate them using a linear reformulation. I created binary variables $y_{p,n,s,i}$ to replace the above variables $c$ and added the three following constraints: $$\begin{align} y_{p,n,s,i} &\leq x_{p,s,i-1}, \forall (p, n) \in C, \forall s \in S, \forall i\in \{1,2,3\} \\ y_{p,n,s,i} &\leq x_{n,s,i}, \forall (p, n) \in C, \forall s \in S, \forall i\in \{1,2,3\} \\ y_{p,n,s,i} &\geq x_{p,s,i-1} + x_{n,s,i} - 1, \forall (p, n) \in C, \forall s \in S, \forall i\in \{1,2,3\} \\ \end{align}$$
However, after doing so, the problem is even slower to solve on the same large instances mentioned before. The issue is that there are many more variables now (one order of magnitude more) and many more constraints. Gurobi, which is the solver I'm using, needs to spend a lot of time in the presolve phase making the problem smaller. Then, when solving the problem, it also spends a large amount of time in the root node and after that starts exploring the tree, where there are now many more nodes to explore, and thus it takes more time to reach the optimal solution.
I tried using a lot of different parameters of Gurobi to no avail. I also added my own cuts with information on the $y$ variables and on the model to help Gurobi but again I didn't succeed in speeding up the optimization.
Is there any other reformulation I could use to remove the bilinear constraint?