The idea given by @Erwin_Kalvelagen is correct, however, if you really want to keep constraints linear we need to avoid product of two different variables appear in our formula. Here is a way to avoid using the deltas. Let's say I have two variables x and y and I want the constraint x≠y. So we need x<y or x>y. x>y is same as x-y=M for a positive M, so we add x-y-M=0 if your linear solver accepts equality for constraints (for example linprog from scipy.optimize) or x-y-M>=0 and putting the bound M>=1. But that's not the end, because if it was, then why not simply adding x-y>=0 without need of an extra auxiliary variable, M, right? Then for x<y, we use N, another positive number. x<y is same as x-y=-N or x-y+N<=0 together with N>=1. But we don't want both of them be satisfied together, because they can't. So what to do? I modify the bounds for M and N to be free (no lower or upper bound) instead of positive, but adding a new constraint M+N>=1, this will guarantee that at least one of them is positive, and because we know both can't be positive at the same time, it means exactly only one will become positive. Thus either x<y or y>x will happen and we have x≠y for sure.

To wrap up x≠y is the same as

x-y-M>=1,
x-y+N<=-1,
M+N>=1.

Of course you can play with it for example using y-x-N>=1 instead of the second line etc. Also I considered integer domain, you can play with it for the other domain scenarios.