Convert system of linear equations to LCP

Question

Background

I'm working on a problem based in physics (constrained dynamics) where I currently have a system of linear equations. I would like to be able to add a few inequality constraints (e.g. for contact constraints), so my thought was to convert the existing system to an LCP and then add constraints there. This also seemed reasonable because I've seen some discussion about Sequential Impulses (a common game physics solver) being equivalent to solving an LCP (e.g. this discussion). However, if this already seems incorrect or that there's a better way to do this, I'm definitely open to suggestions.

I have tried the following steps:

1. Observe that $A_0 + A_1x_1 + \cdots + A_nx_n=0$ is equivalent to $$A_0 + A_1x_1 + \cdots + A_nx_n \ge 0 \land A_0 + A_1x_1 + \cdots + A_nx_n \le 0$$ which is itself equivalent to $$A_0 + A_1x_1 + \cdots + A_nx_n \ge 0 \land -A_0 - A_1x_1 - \cdots - A_nx_n \ge 0$$
2. Observe that $x_i \in \mathbb{R}$ can be substituted for $u_i - v_i$ where $u_i, v_i \ge 0$.
3. Starting from $Ax - b = 0$ then, we can form an LCP where $$M = \begin{bmatrix} A_{11} & -A_{11} & \cdots & A_{1n} & -A_{1n} \\ -A_{11} & A_{11} & \cdots & -A_{1n} & A_{1n} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ A_{n1} & -A_{n1} & \cdots & A_{nn} & -A_{nn} \\ -A_{n1} & A_{n1} & \cdots & -A_{nn} & A_{nn} \end{bmatrix}\qquad q = \begin{bmatrix}-b_1\\ b_1\\ \vdots\\ -b_n\\ b_n \end{bmatrix} $$ The solution to the original system becomes $$x = [z_1 - z_2, \ldots, z_{n-1} - z_n]^T$$

Example

Consider the intersection of three planes: $$x - y = 0\\ y + z = 2\\ x - z = 0$$ As a linear system, this is $$A = \begin{bmatrix} 1 & -1 & 0 \\ 0 & 1 & 1 \\ 1 & 0 & -1 \end{bmatrix}\qquad b = \begin{bmatrix}0 \\ 2 \\ 0 \end{bmatrix}$$ The solution is $x = [1, 1, 1]^T$.

Using my proposed LCP transformation, we get $$M = \begin{bmatrix} 1 & -1 & -1 & 1 & 0 & 0 \\ -1 & 1 & 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & -1 & 1 & -1 \\ 0 & 0 & -1 & 1 & -1 & 1 \\ 1 & -1 & 0 & 0 & -1 & 1 \\ -1 & 1 & 0 & 0 & 1 & -1 \end{bmatrix}\qquad q = \begin{bmatrix}0\\ 0\\ -2\\ 2\\ 0\\ 0\end{bmatrix} $$ One solution to this is $z = [1, 0, 1, 0, 1, 0]^T \quad w = 0$, which produces the correct $x$ for the original system.

Problem

I have tested the above example with Siconos and lemkelcp. Both fail to converge on a solution, with the latter specifically returning "secondary ray found". I don't know what to make of this.

Question

So my questions are:

1. Is there a better or commonly known transformation from a system of linear equations to an LCP?

2. What might be a reason that numerical algorithms for LCPs fail to converge for this transformation and can that be fixed? The two things that are suspicious to me are that $x = u - v$ introduces infinite solutions and I'm not sure if my neglecting to consider $w$ causes a bad situation.

3. As mentioned earlier, is there another approach I might take to introduce inequality constraints into a linear system?

I'm pretty out of my depth here mathematically, so any help is greatly appreciated. Thanks!