For linear programming type bounds it is sometimes only possible to give effective bounds (that is bounds that work and are manageable) and what is surprising is that the primitive method often gives in fact optimal bounds.
The LP problem from McEliece, Rodemich, Rumsey and Welch's paper that is cited in the question requires the auxiliary function $\beta(x)$ satisfy $\beta(j) = 1+ \sum_{k=1}^{n} y_k K_k (j)\leq 0$ for all $j=d,d+1,\dots,n$. The supplied $\beta(x)$ is designed to meet this requirement by making it change sign at value $a$, so that $\beta(x)\leq 0$ for $x\geq a$. This is the first point.
The Krawtchouk polynomials already appear in the LP problem, so that there are no questions of why the might appear in the auxiliary function $\beta$, but just to emphasize the importance of the Krawtchouk polynomials, they are used in discrete linear programming problems due to the positive definiteness criterion associated to them, namely that for a polynomial $$f(z)=\sum\limits_{i=0}^{n} a_{i}z^{i}$$ the matrix $$f(d(x,y)),\ x,y\in\mathbb{F}^{n}$$ is non-negative definite if and only if all coefficients $\lambda_{i}$ of the expansion $f(z)=\sum\limits_{i=0}^{n} \lambda_{i} K_{i}(z)$ over Krawtouck polynomials are nonnegative.
Now, to keep all the coefficients $y_k$ positive as required in the LP program, it makes sense to introduce the square, but in order to preserve the sign change in $\beta(x)$ at $x=a$, one divides by $(a-x)$.
Finally, the choice of the exact expression inside the square works, because as was noted, the Christoffel-Darboux formula allows for rewriting
$$K_t(a) K_{t+1}(x) - K_{t+1}(a)K_{t}(x)=\frac{2(a-x)}{t+1}\binom{n}{t}\sum\limits_{k=0}^t \frac{K_{k}(x)K_{k}(a)}{\binom{n}{k}}$$
so that one may check quickly that $\beta(x)$ has expansion coefficients $y_{k}$ in the Krawtouck polynomials non-negative. Optimizing in $a$ and $t$ as noted give the MRRW upper bound $M_{LP}(n,d)\leq \binom{n}{t}\frac{(n+1)^2}{2(t+1)}$.