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, ...

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,\dots,y_n]\cdot A^l\cdot [1,1,\dots,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, ...