How about reformulating your problem as QUBO by including the constraint into the objective function and treating $x_i$ as binary and then doing some post-processing with measured samples. The constraint of your problem allows to exploits some tricks.
The state $| \psi \rangle$ representing an optimized QAOA circuit will have the following form \begin{align} | \psi \rangle = a_0 |0 \rangle + \cdots +a_{2^n-1} | 2^{n}-1 \rangle. \end{align} Given that the circuit is optimized, it is expected that the state $|k \rangle$ has a probability $|a_k|^2$ near zero for any $k \notin \{1, 2, 4,...,2^{n-1}\}$. This is because for $k \notin \{1, 2, 4,...,2^{n-1}\}$ we get states like $|3 \rangle = |0 \cdots011\rangle$ or $|5 \rangle = |0 \cdots101\rangle$. Such states violate the constraint and introduce penalty. Since QAOA tries to minimize the expected energy it will try to decrease the probability $|a_k|^2$ of such infeasible states. Therefore, the most probable states are $|k\rangle$ for $k \in \{1, 2, 4,...,2^{n-1}\}$. Or in other words, states that have only single bit on e.g. $|1\rangle = |0\cdots01\rangle$ or $|2\rangle = |0\cdots10\rangle$.
Next, you can make $m$ measurements of the circuit to obtain samples. For simplicity, suppose that your problem has 3 variables. Then you sample set will contain many binary strings of the form $001, 010, 100$. These are all feasible solutions. Of course there will be infeasible strings too. But you can ignore them.
Finally, select all the feasible strings and for each bit compute the frequency of 1. This will give the "averaged" versionestimate of a solution which always satisfies the constraintprobability distribution given by $|\psi\rangle$ for each feasible solution. For example, suppose you sampled 5 bit strings $100, 100, 100, 010, 001$ then the frequency vector for each bit is $(3/5, 1/5, 1/5)$. This is the approximate solution to the problemapproximation of probabilities given by $| \psi \rangle$ that only includes feasible solutions. And we always get $x_1 + x_2 + x_3 = 1$.
In the other post it was suggested to use a binary expansion for $x_i$ so that it could model float numbers with limited precision. While it is something that people usually do, such an approach will make your problem extremely hard in comparison to the original problem. First, with each variable $y_k$ you double the search space (and increase the number of qubits). Second, the energy landscape becomes really rugged. Basically, you turn an "easy" problem into a really hard one.