In purely technical terms, when we have a perturbation problem we introduce a dimensionless parameter $\lambda$ and rewrite the problem as unperturbed part $+\lambda$perturbation. Then we expand the solution in power series of $\lambda$ and equate coefficients. Using your concrete example to demonstrate, we have
$$y''-f(x)y'-w(x)y/x=0$$
as the unperturbed equation, to which we know a set of solutions (otherwise it is useless to use as an unperturbed problem) $y_n(x)$ and let's focus on one solution $y_0(x)$. Now for the perturbation we now expand $Y(x) = y_0(x) + \sum_{n>0}\lambda^n q_n(x)$ and we plug this into the perturbed equation
$$Y''-f(x)Y'-w(x)Y/x+\lambda\epsilon[g(x)-w(x)Y/x]=0$$
and we equate powers of $\lambda$. Then at the end we set $\lambda=1$ and we get the solution. Throughout all this process, we didn't assume anything on the perturbation, besides the fact that it can be expanded in a power series (which is not trivial! A lot of functions do not have such an expansion.) We also didn't have any dimensional discussion - $\lambda$ is dimensionless. However, as you can see in your example, $\lambda$ always came with $\epsilon$, so the power series in $\lambda$ will also be a power series in $\epsilon$ (this is not always the case, where we can in advance identify a parameter to expand in).
Now, we must analyze our solution, and see what we get. The solution might be a complete nonsense, of course, and no one guarantees that our approximation method works. In order for the solution to be valid, we want that the power series will converge. So we look at the results $q_n(x)$ and want them to be smaller than the original result itself $||q_n(x)|| \ll ||y_0(x)||$ (note that as $\lambda$ is dimensionless these two properties are bound to have the same dimensions). But it might be that this is not necessarily the "smallness" we are looking for. It might be that we want our solutions to remain well separated from the other solutions, in some sense $||q_n|| \ll ||y_0-y_n||$. Or we might want some observable property to get a small relative correction $| E[y_0+\sum_{n>0} q_n]-E[y_0] | \ll E[y_0]$.
If we have a small dimensionful expansion parameter in advance, as in your case with $\epsilon$, then it should correspond to some scale in our solution $y_0$, and then we can compare them and assume that for the most cases if $|\epsilon_{y_0}| \gg |\epsilon|$ then the expansion in power-series of $\epsilon$ will represent small corrections to the true solution. We can make it more formal by dividing the equation with $\epsilon_{y_0}$ and get a power-series in a dimensionless small number $\epsilon/\epsilon_{y_0}$. This, however, depends on the particulars of the problem.
By the way, it might be that the power-series does not converge, and still we will use the first terms in it as a good approximation. This is called an asymptotic expansion.
How well this approximation works depends on how accurate we want our result to be. In general if we expand to $(\epsilon/\epsilon_{y_0})^n$ then this is the level of accuracy of our result. Sometimes just the first or two leading orders are enough for what we want, sometimes not.
