Farkas' Lemma is indeed the way to go, but we need the right setting. Below I give a sketch.
For simplicity, assume that we work at a vertex $x=0$ of $P$.
So we want to find a minimal set of generators for the cone $\DeclareMathOperator{\cone}{cone}C:=\cone(P)=\cone (\mathcal V)$, where $\mathcal V\subseteq P$ is the set of vertices of $P$.
What we want to understand is whether every such "minimal generator" $y\in\mathcal V$ is a neighbor of $x$, because if so, then the edge-directions indeed generate $C$.
So, suppose that $y\in \mathcal V$ is part of such a minimal set of generators. Then $y\not\in C':=\cone(\mathcal V\setminus \{y\})$ (here you need to use that no three vertices of $P$ are colinear). By Farkas' Lemma, we can then separate $y$ from $C'$ via a hyperplane. In particular, we can choose this hyperplane with normal vector $n$ so that
$$\def\<{\langle}\def\>{\rangle}\<n,x\>=0,\quad\<n,y\> >0\quad\text{and}\quad\<n,z\><0\text{ for all $z\in \mathcal V\setminus\{x,y\}$}.$$
It is not too hard to argue that we can choose $n$ linearly independent from $y$ (if we are working in dimension $d\ge 2$). Then
$$n':=n-y\frac{\<n,y\>}{\<y,y\>} \not=0.$$
You can check that we have $\<n',x\>=\<n',y\>=0$ and $\<n',z\><0$ for all $z\in \mathcal V\setminus\{x,y\}$ (the latter needs some thought, but is possible). In other words, the hyperplane orthogonal to $n'$ supports $P$ exacty at the two vertices $x$ and $y$, which proves that these form an edge of $P$.
In still other words, $\cone(P)$ is generated by the neighbors of $x$.
Some further explanation
As requested in the comments, I elaborate on $\<n',z\><0$ for all $z\in\mathcal V\setminus\{x,y\}$. As Epiousios noted, this is the same as
$$(*)\quad \underbrace{\<n,z\>}_{<0} < \underbrace{\frac{\<n,y\>}{\<y,y\>}}_{>0} \<y,z\>,$$
which would be obviously true if $\<y,z\>>0$. However, this is not always the case.
But, we can do a trick: before we start with any of our argument, we can tranform our polytope $P$ into an more convenient polytope $P'$, for which any two neighbors $y,z$ of $x=0$ satisfy $\<y,z\>>0$ (meaning $\sphericalangle(y,z)<90^\circ$).
We can do this by stretching $P$ in a certain way. Hopefully, the following image makes this clearer:
Since this is a linear transformation, this changes nothing about the actual problem. But this time $(*)$ is trivially satified.
Let's assume wlog. $P$ has full rank by restricting ourselves to the affine hull if necessary. Let $x$ be a point in the interior of $P$.
Claim: For every $i$, there exists a convex combination $x = \sum_j \lambda^i_j x_j$ such that $\lambda^i_i \neq 0$.
Then $x = \sum_j \frac 1n(\sum_i\lambda^i_j)x_j$ is a positive convex combination of $x$, as desired.
Proof of claim: Since $x$ is in the interior of $P$, there exists some $\varepsilon > 0$ such that $y := x + \varepsilon(x - x_i)$ is in $P$. Hence there exists a convex combination $y = \sum_j\mu_jx_j$. Then $$x = \frac{1}{1+\varepsilon}(y + \varepsilon x_i) = \frac{\varepsilon}{1+\varepsilon}x_i + \sum_j\frac{1}{1+\varepsilon}\mu_j x_j,$$ so
$$\lambda_j^i := \frac{1}{1+\varepsilon}\begin{cases}\varepsilon+\mu_j&\text{ if }i=j\\\mu_j&\text{ else }\end{cases}$$
satisfies $x = \sum_j \lambda^i_j x_j$ and $\lambda^i_i \neq 0$.
Best Answer
You do not need to express the convex hull explicitly as a system of inequalities $\mathbf{A x} \leq \mathbf{b}$. The convex hull, by definition, is $$ \operatorname{conv}(\mathbf{v}_1, \dots, \mathbf{v}_m) = \left \{ \sum_{i=1}^m w_i \mathbf{v}_i : \sum_{i=1}^m w_i=1, w_i \geq 0 \right\} $$
Therefore, the projection onto the convex hull of a vector $\mathbf{y}$ can be computed from the solution of $$ \min_{\mathbf{w}} \quad \left \| \sum_{i=1}^m w_i \mathbf{v}_i - \mathbf{y} \right \|_2^2 \quad \text{s.t.} \quad \sum_{i=1}^m w_i=1, w_i \geq 0 $$ The vector $\mathbf{w}$ contains the barycentric coordinates of the projection. I am not aware of a specialized algorithm of computing the optimal $\mathbf{w}$, but any QP solver can do so.
A simple solution method which you can implement yourself can be obtained from a fast projected gradient algorithm ,such as FISTA, and remembering that the constraints are exactly the unit simplex, and there is a simple $O(m \log m)$ algorithm for projecting a point onto the unit simplex. Moreover, you also have the well-known Mirror Descent algorithm for optimizing over the unit simplex, which can be quite fast in practice.