I am trying to solve a problem that involves constraints in which products of two decision variables appear. So far, I read that such products can be reformulated to a difference of two quadratic terms:
$$ x_1\cdot x_2=y_1^2-y_2^2 $$
where $y_1 = 0.5 \cdot \left( x_1+x_2 \right)$ and $y_2 = 0.5 \cdot \left( x_1-x_2 \right)$. As stated in "Model building in mathematical programming" by H. P. Williams, I tried to linearize $y_1^2$ and $y_2^2$ by piecewise approximation. When I only use two node points (one interval) for the piecewise approximation, my solvers (Gurobi 5.6.3 and CPLEX 12.5.1) are able to solve the problem, but when I introduce more node points, both conclude that the problem becomes infeasible.
I already tried SOS2 variables as well as a binary approach for the linearization of $y_1^2$ and $y_2^2$.
I also used Taylor series to directly approximate $x_1\cdot x_2$ which resulted in similar results as the approximation using two node points. Since this method succeeded, I suppose that the rest of my model is feasible and only this product makes things complicated.
What can cause the infeasibility? I thought that adding more node points would be beneficial to the problem.
Are there better ways to approximate such products? I have also experimented with a reformulation using logarithms, but this caused more computational effort and did not solve the infeasibility problem.
The AIMMS Modelling Guide in section 7.7 provides the following hint for bounded variables: