As you suggested, your problem can be formulated as a quadratic program:
$$\boxed{\begin{array}{rl} \min & \|\mathbf x_0-(\alpha_1\mathbf v_1+\cdots +\alpha_n\mathbf v_n)\|^2\\ \text{s.t.} & \alpha_1+\cdots+\alpha_n=1\\ & \alpha\ge 0 \end{array}}$$
The optimal solution to this program gives you the projection $\mathbf v:=\alpha_1 \mathbf v_1+\cdots \alpha_n\mathbf v_n$ in barycentric coordinates.
You can express above program in a more "standard" way, by setting $V := (\mathbf v_1,...,\mathbf v_n)\in\Bbb R^{d\times n}$ to be the matrix with the $\mathbf v_i$ as columns, and optimize $\mathbf v\in\Bbb R^d,\alpha\in\Bbb R^n$ via
$$\boxed{\begin{array}{rl} \min & \|\mathbf x_0-\mathbf v\|^2\\ \text{s.t.} & V\alpha = \mathbf v\\ & \alpha\ge 0 \end{array}}$$
In the optimal point, $\mathbf v$ is the desired projection and $\alpha$ contains the barycentric coordintes. There exist standard solvers for quandratic problems like this.