Coming a bit late to the party, but since there is no answer, let me offer one.

Yes there is a program to make progress towards the YM problem, although with no guarantee of success in the end. The methods/part of mathematics that have the potential to lead to working this program out, that's called constructive quantum field theory (CQFT). To learn about it, the best (although not perfect) reference that I can recommend is the book <a href="http://www.rivasseau.net/resources/book.pdf">"From Perturbative to Constructive Renormalization"</a> by Vincent Rivasseau. A brief account is <a href="https://arxiv.org/abs/1102.4238">"Mode d'emploi de la th&eacute;orie constructive des champs bosoniques (A user's guide to bosonic constructive field theory)"</a> by J&eacute;r&eacute;mie Unterberger, but it is in French.

**The program of study:**

Solve, in the order (of difficulty) stated below, the analogue of the YM construction and mass gap problem for the following two models.

 1. The two-dimensional Gross-Neveu model with large (but finite) number $N$ of components.
 2. The two-dimensional sigma-model, also with large $N$.

These two problems are within reach of existing CQFT methods. There is a very steep increase in difficulty from these two models to 4d YM, so there is no guarantee that if one solves these problems one would be able to also do YM. However, in my opinion, there is no way one can solve YM if one cannot solve the above two problems.

The YM Clay Millenium problem requires constructing the UV limit in infinite volume and showing a mass gap. There is a lighter version of that, which also essentially is a necessary milestone along the way. It is to look at the lattice problem 
without tackling the UV limit. Namely, the question is to rigorously control the infinite volume limit on the unit lattice but at *arbitrarily low temperature*.
One often reads in the physics literature that there is no problem doing things rigorously on the lattice. This is false if one is talking about the infinite lattice. At high temperature, this is nontrivial but can be done using cluster expansions.
For YM and the sigma-model (Fermions also behave like compact-valued random variables) the random variables take values in compact spaces, so by a variant of Prokhorov's Theorem, one can construct subsequential infinite volume limits. However, these are almost useless. If one wants to study this infinite volume problem on the lattice, a top priority intermediate question that should be settled first is that of establishing ferromagnetic correlation inequalities. There is evidence that non-Abelian models like the $O(N)$ model (the lattice version of the sigma-model) or lattice gauge theories do satisfy correlation inequalities of ferromagnetic type.