Is there a way to use Euler-Lagrange equations to derive the following formula?
I can derive the term not due to damping, I guess the damping force should be somehow modeled using a generalized force. Can you please show me how?
My attempt was starting by
$$ F_{d,ij} = k_d \lVert v_{ij} \rVert u_{ij} $$
Where $v_{ij}$ is the relative velocity vector and $u_{ij}$ is the unit vector going from mass $i$ to $j$. But I'm missing how should I go from here to the generalized forces.
Thank you.

OP is essentially asking to find a velocity-dependent potential for coupled damped oscillators. This can be reduced to the question of finding a velocity-dependent potential $U({\bf x},\dot{\bf x},t)$ for the 1-particle force $$ {\bf f} ~=~-k_d{\bf x}\frac{{\bf x} \cdot\dot{\bf x}}{|{\bf x}|^2} ~=~-k_d{\bf x}\frac{d\ln |{\bf x}|}{dt}.$$ One may show, using methods e.g. outlined in my Phys.SE answer here that such velocity-dependent potential $U({\bf x},\dot{\bf x},t)$ does not exist. So OP's model cannot be obtained from a stationary action principle (under the assumption that we shouldn't modify the kinetic sector of the action).
However, it is still possible to describe OP's model via Lagrange equations $$ \frac{d}{dt} \left(\frac{\partial L}{\partial \dot{\bf x}_i}\right)-\frac{\partial L}{\partial {\bf x}_i}~=~-\frac{\partial {\cal F}}{\partial \dot{\bf x}_i} $$ using the Rayleigh dissipation function$^1$ $$ {\cal F} ~=~ \frac{1}{2}\sum_{i<j} k_d^{ij} \frac{({\bf x}_{ij} \cdot\dot{\bf x}_{ij})^2}{|{\bf x}_{ij}|^2}, \qquad {\bf x}_{ij}~:=~{\bf x}_{i}-{\bf x}_{j}. $$ The Lagrangian itself is $$L~=~T-V, \qquad T~=~\frac{1}{2}\sum_i m_i |\dot{\bf x}_i|^2, \qquad V~=~\frac{1}{2}\sum_{i<j}k^{ij}_s (|{\bf x}_{ij}| - \ell^{ij}_0)^2. $$ See also this related Phys.SE post.
--
$^1$ In this answer, we have slightly generalized OP's model to allow the spring and dissipation constants to depend on the particle pair $(i,j)$.