Actually, at least for a single point subjected to a friction force $F= -\gamma v$ and other forces associated with a potential $U(t,x)$ there exists a Lagrangian: $${\cal L}=e^{t\gamma/m}\left(\frac{m}{2}\dot{x}^2 -U(t,x)\right)\:.$$
The point is that this Lagrangian is not of the form $T-U$, nevertheless it gives rise to the correct equation of motion, the same obtained by using the Rayleigh dissipation function you mentioned.
This Lagrangian however cannot be produced by direct application of the principle of virtual works you mention.
I) The Euler-Lagrange (EL) equations behave covariantly under reparametrizations$^1$ of the form
$$ \tag{1} q^{\prime i}=f^i(q,t),$$
i.e. it is equivalent to reparametrize before or after forming the EL equations.
II) The above property even holds for a Lagrangian $L(q,\dot{q},\ddot{q},\ldots, \frac{d^Nq}{dt^N};t)$ that depends on higher-order time-derivatives, although a higher-order version of Euler-Lagrange equations with higher-order derivatives is needed in such case.
III) However, for a velocity-dependent reparametrization $q^{\prime }=f(q,\dot q,t)$, which OP mentions in his second line (v2), the substitution before or after in general leads to EL eqs. of different orders. We expect that the higher-order EL eqs. to always factorize via the corresponding lower-order EL eqs., so that solutions to the lower-order EL eqs. are also solutions to the higher-order EL eqs. but not vice-versa.
Similarly for acceleration-dependent reparametrizations, etc.
IV) Example: Consider the velocity-dependent reparametrization
$$\tag{2} q^{\prime}~=~q+A \dot{q}, \qquad A>0,$$
of the Lagrangian$^2$
$$\tag{3} L^{\prime}~=~ \frac{1}{2} q^{\prime 2}~=~\frac{1}{2}(q+A \dot{q})^2~\sim~ \frac{1}{2}q^2 +\frac{A^2}{2} \dot{q}^2. $$
(We call $q^{\prime}$ and $q$ the old and new variables, respectively.)
Before, the EL equation is of first order in the new variables$^3$
$$\tag{4} 0\approx q^{\prime}~=~q+A \dot{q},$$
with only exponentially decaying solutions. After the reparametrization, the EL equation is of second order
$$\tag{5} 0\approx q- A^2 \ddot{q}~=~(1-A\frac{d}{dt})(q+A \dot{q}),$$
so that it has more solutions. Note however that eq. (5) factorize via (=can be obtained from) eq. (4) by applying a differential operator $1-A\frac{d}{dt}$.
--
$^1$ There are various standard regularity conditions on a reparametrization (1) such as e.g. invertibility and sufficiently differentiability. The higher jets (velocity, acceleration, jerk, etc) are implicitly assumed to transform in the natural way.
$^2$ The $\sim$ sign means here equal modulo total derivative terms.
$^3$ The $\approx$ sign means here equal modulo the EL equations.
Best Answer
If the force is not derived from a potential, then the system is said to be polygenic and the Principle of Least Action does not apply. However, the Euler-Lagrange equations can be derived from d'Alembert Principle.
If we decompose the applied (or specified) forces acting on particle $\alpha$ into monogenic (derived from a potential), $\vec F_\alpha^m$ and polygenic forces, $\vec F_\alpha^p$, then d'Alembert Principle reads, $$\sum_\alpha(\vec F_\alpha^m+\vec F_\alpha^p-\dot{\vec p}_\alpha)\cdot\delta\vec r_\alpha=0.$$ The next step is to write this equation in terms of generalized coordinates $q_i$. The result is the following equation of motion $$\frac{d}{dt}\frac{\partial T}{\partial \dot q_i}-\frac{\partial T}{\partial q_i}=Q_i^m+Q_i^p,$$ where $$Q_i^p\equiv\sum_\alpha\vec F_\alpha\cdot\frac{\partial \vec r_\alpha}{\partial q_i}.\tag1$$
The monogenic force can be obtained from a potential $V$, $$Q_i^m=-\frac{\partial V}{\partial q_i},$$ hence the equation of motion $$\frac{d}{dt}\frac{\partial T}{\partial \dot q_i}-\frac{\partial T}{\partial q_i}+\frac{\partial V}{\partial q_i}=Q_i^p.$$ If the potential does not depend on velocities, then this equation can also be written as $$\frac{d}{dt}\frac{\partial L}{\partial \dot q_i}-\frac{\partial L}{\partial q_i}=Q_i^p,\tag2$$ where $L=T-V$ is the Lagrange function. Equation (2) is the one you shall use, together with Eqn. (1) to obtain the generalized force $Q_i^p$.
Edit:
Let's now apply this approach to the example posed in the question. There are two external forces, which can be written as $\vec F_1 = [F_x(t) \; , 0]^T$ and $\vec F_2 = [0 \; , -F_\theta(t)]^T$. The position of each body (regarded as a point mass) is $\vec r_1 = [x \; , 0]^T$ and $\vec r_2 = [x + l \sin \theta \; , -l \cos \theta]^T$. Therefore, we calculate $$Q_1^p = \vec F_1 \cdot\frac{\partial \vec r_1}{\partial x} + \vec F_2 \cdot\frac{\partial \vec r_2}{\partial x} = F_x(t)$$ and $$Q_2^p = \vec F_1 \cdot\frac{\partial \vec r_1}{\partial \theta} + \vec F_2 \cdot\frac{\partial \vec r_2}{\partial \theta} = -F_\theta(t) l \sin \theta .$$
Finally, the corresponding Euler-Lagrange equations are $$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot x}\right )-\frac{\partial L}{\partial x}= F_x(t)$$ $$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot \theta}\right )-\frac{\partial L}{\partial \theta}=-F_\theta(t) l \sin \theta,$$ where $$L = T - V = \frac{M}{2} \left\lVert\dot {\vec r_1}\right\rVert^2 + \frac{m}{2} \left\lVert\dot {\vec r_2}\right\rVert^2 + m g l \cos \theta .$$