Sisco,
You ask about modeling a 0/1 variable using smooth constraints.
If all you want is x to be binary, you could do x * (1-x) = 0. This won't be convex but that's to be expected if you have binary variables.
There is some relationship between 0/1 variables and complementarity, but they are not the same thing. Usually a model can be more directly put into one form or the other. For example, a model with a lot of yes/no decisions (e.g. do we build a warehouse here or not?) is naturally modeled with binary variables. A model with regime changes or phase changes (e.g. production is set to satisfy an optimality condition, as long as it is above the contractual or physical lower bound) is more naturally done with complementarity. Also there is the whole weight of previous work - if you use binary variables in cases where most others do, you'll likely have a wealth of examples and help available. Ditto for complementarity.
If you throw in the possibility of approximating the complementarity with a sigmoid or other smooth function (there are several options in GAMS, see the table here
https://www.gams.com/latest/docs/UG_Par ... ch=sigmoid), you introduce other issues. The problem may get easier to solve and you can use different solvers, but you solve only an approximation. As the approximation gets better, it also gets more difficult: it will have a tighter kink and perhaps other scaling issues.
It isn't really a matter of solver preference as much as capability. If you use binary variables, you must use a solver that can handle discrete variables. If you use complementarity, you have to use an MCP or MPEC solver. If you use smoothing functions, you can use an NLP solver. There's not a lot of overlap between these three.
HTH,
-Steve
