This is quite simple: if $\|f - g\|_\infty \leqslant \epsilon$, then clearly
$$ \sum_{i = 1}^n |f(x_i) - f(x_{i-1})| \leqslant \sum_{i = 1}^n |g(x_i) - g(x_{i-1})| + 2 n \epsilon \leqslant V(g) + 2 n \epsilon.$$
Therefore,
$$ \sum_{i = 1}^n |f(x_i) - f(x_{i-1})| \leqslant \inf_{g \in C_\epsilon(f)} V(g) + 2 n \epsilon,$$
and consequently
$$ \sum_{i = 1}^n |f(x_i) - f(x_{i-1})| \leqslant \liminf_{\epsilon \to 0^+} \inf_{g \in C_\epsilon(f)} V(g) .$$
Take the supremum over all possible partitions to get
$$ V(f) \leqslant \liminf_{\epsilon \to 0^+} \inf_{g \in C_\epsilon(f)} V(g) .$$
EDIT: I forgot about the converse inequality:
On the other hand, $f$ can be approximated by functions from $C_\epsilon(f)$ in the total variation norm (see below), so that
$$ V(f) \geqslant \inf_{g \in C_\epsilon(f)} V(g) , $$
and hence
$$ V(f) \geqslant \limsup_{\epsilon \to 0^+} \inf_{g \in C_\epsilon(f)} V(g) . $$
Regarding approximation: if $f$ is non-decreasing, then there is a smooth function $g \in C_\epsilon(f)$ such that $V(g) \leqslant V(f)$ (just mollify $f$ appropriately). In the general case, if $V(f) < \infty$, then $f = f_+ - f_-$ with $f_+$ and $f_-$ non-decreasing and $V(f) = V(f_+) + V(f_-)$. Find $g_+$ and $g_-$ as above; then $g = g_+ - g_-$ is in $C_{2\epsilon}(f)$, and $$V(g) \leqslant V(g_+) + V(g_-) \leqslant V(f_+) + V(f_-) = V(f).$$