I'm not sure if the formula below is useful, but at least it can be used to numerically count edge-invariant walks classes for small $n,l$.

Let's assign a unique variable to each edge and each vertex of the complete graph, say, $x_{ij}$ to an edge $(i,j)$ and $y_i$ to a vertex $i$. Then, we assign to each edge $(i,j)$ the weight $x_{i,j}y_j$ and consider the graph adjacency matrix $A$. For example, for $n=4$ the matrix is
$$A = \begin{bmatrix} 
x_{11}y_1 & x_{12}y_2 & x_{13}y_3 & x_{14}y_4\\
x_{21}y_1 & x_{22}y_2 & x_{23}y_3 & x_{24}y_4\\
x_{31}y_1 & x_{32}y_2 & x_{33}y_3 & x_{34}y_4\\
x_{41}y_1 & x_{42}y_2 & x_{43}y_3 & x_{44}y_4
\end{bmatrix}.$$

Then the number of classes of edge-invariant of length $l$ is given by the number of distinct monomial terms in
$$[y_1,y_2,y_3,y_4]\cdot A^l\cdot [1,1,1,1]^T.$$

For example, for $n=2$ and $l=3$ we get the polynomial:
$$x_{11}^3y_1^4 + x_{11}^2x_{12}y_1^3y_2 + x_{11}^2x_{21}y_1^3y_2 + 2x_{11}x_{12}x_{21}y_1^3y_2 + x_{11}x_{12}x_{21}y_1^2y_2^2 + x_{12}^2x_{21}y_1^2y_2^2 + x_{12}x_{21}^2y_1^2y_2^2 + x_{11}x_{12}x_{22}y_1^2y_2^2 + x_{11}x_{21}x_{22}y_1^2y_2^2 + x_{12}x_{21}x_{22}y_1^2y_2^2 + 2x_{12}x_{21}x_{22}y_1y_2^3 + x_{12}x_{22}^2y_1y_2^3 + x_{21}x_{22}^2y_1y_2^3 + x_{22}^3y_2^4,$$
where the monomials describe equivalence classes of walks (more specifically, $x$'s describe $e_p$ and $y$'s describe $f_p$) and the coefficients give the size of each class. There are $14$ distinct monomials here and only two of them have coefficient $2$, as expected.

Here is my SAGE implementation of this formula:


    # generator for variables
    class VariableGenerator(object): 
      def __init__(self, prefix): 
         self.__prefix = prefix 
      @cached_method 
      def __getitem__(self, key): 
         return SR.var("%s%s"%(self.__prefix,key)) 

    def NumEIClasses2(n,l):
      x = VariableGenerator('x')
      y = VariableGenerator('y')

      R = PolynomialRing(QQ,[x[i] for i in range(n*n)] + [y[i] for i in range(n)])

      A = matrix([[x[n*i + j]*y[j] for j in range(n)] for i in range(n)])

      u = vector([1 for i in range(n)])
      Y = vector([y[i] for i in range(n)])

      P = R( (Y.row() * A^l * u.column())[0,0] )
      #print P    # print the resulting polynomial

      return P.hamming_weight()

For example, for $n=4$ it gives counts 16, 64, 244, 856, 2728, 7892, 20876, 51020, 116408, ...