# Failure to solve system of two exponential equations in 2 unknowns

ClearAll[xaa, xab, xba, xbb]
A = 6349;
B = -384;
R = .08314;
T = 323;
P = 1;
NSolve[{xaa*Exp[(A + 3*B)/(R*T)*xba^2 - 4*B/(R*T)*xba^3] -
xab*Exp[(A + 3*B)/(R*T)*xbb^2 - 4*B/(R*T)*xbb^3] == 0,
xba*Exp[(A - 3*B)/(R*T)*xaa^2 + 4*B/(R*T)*xaa^3] -
xbb*Exp[(A - 3*B)/(R*T)*xab^2 + 4*B/(R*T)*xab^3] == 0,
xaa + xba == 1,
xab + xbb == 1},
{xaa, xab, xba, xbb}, Reals]


I am trying to solve a liquid liquid equilibrium problem from thermodynamics. But when I use this code I get the following as output:

NSolve[{E^(193.526 xba^2 + 57.1977 xba^3) xaa -
E^(193.526 xbb^2 + 57.1977 xbb^3) xab == 0,
E^(279.323 xaa^2 - 57.1977 xaa^3) xba -
E^(279.323 xab^2 - 57.1977 xab^3) xbb == 0, xaa + xba == 1,
xab + xbb == 1}, {xaa, xab, xba, xbb}, Reals]


If I use Solve instead, it complains about me trying to solve a system with inexact coefficients. Some other constraints I have added to the list of functions with no luck include: 0<=xaa<=1 0<=xab<=1 0<=xba<=1 0<=xbb<=1

ClearAll[xaa, xab, xba, xbb]
A = 6349;
B = -384;
R = .08314;
T = 323;
P = 1;

eqns = {xaa*Exp[(A + 3*B)/(R*T)*xba^2 - 4*B/(R*T)*xba^3] -
xab*Exp[(A + 3*B)/(R*T)*xbb^2 - 4*B/(R*T)*xbb^3] == 0,
xba*Exp[(A - 3*B)/(R*T)*xaa^2 + 4*B/(R*T)*xaa^3] -
xbb*Exp[(A - 3*B)/(R*T)*xab^2 + 4*B/(R*T)*xab^3] == 0, xaa + xba == 1,
xab + xbb == 1};

(sol = FindRoot[
eqns, {{xaa, #}, {xab, #}, {xba, 1 - #}, {xbb, 1 - #}}] & /@
Range[0, 1, .1]) // Column And @@ (And @@@ (eqns /. sol))

(* True *)


So the general solution is {xaa == xab, xba == xbb == 1 - xaa}

And @@ (eqns /. {xab -> xaa, xbb -> xba} /. xba -> 1 - xaa)

(* True *)