Dado un referencial con vectores base $(\hat i, \hat j,\hat k)$ y realizando
$$ \cases{ q=r(\hat i \cos(\omega t) + \hat j\sin(\omega t))+q_0 \hat k\\ p = (x(t),y(t),z(t)) }\ \ \ \ \ (1) $$
el lagrangiano es
$$ L(p,\lambda, t) = \frac 12m\|\dot p\|^2-m g p\cdot\hat k+\lambda(\|p-q\|^2-l^2) $$
Aquí $\lambda$ es un multiplicador de lagrange para considerar en la variación también la restricción $\|p-q\|=l^2$ y las ecuaciones de movimiento se obtienen de
$$ \cases{ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot p}\right)-\frac{\partial L}{\partial p}=0\\ \frac{d^2}{dt^2}(\|p-q\|^2-l^2) = 0\\ } $$
o
$$ \cases{2 r \lambda \cos (\omega t)-2 \lambda x(t)+m x''(t)=0\\ 2 r\lambda \sin (\omega t)-2 \lambda y(t)+m y''(t)=0\\ g m+2 \lambda q_0-2 \lambda z(t)+m z''(t)=0\\ r \cos (\omega t)\left(\omega ^2 x(t)-2 \omega y'(t)-x''(t)\right)+r \sin(\omega t) \left(\omega ^2 y(t)+2 \omega x'(t)-y''(t)\right)+(z(t)-q_0) z''(t)+x(t)x''(t)+x'(t)^2+y(t) y''(t)+y'(t)^2+z'(t)^2=0 } $$
Resolviendo para $\{x''(t),y''(y),z''(t),\lambda\}$ tenemos
$$ \cases{ x''(t)= \frac{(r \cos (\omega t)-x(t)) \left(r \omega \sin (\omega t) \left(\omega y(t)+2 x'(t)\right)+r \omega \cos (\omega t) \left(\omega x(t)-2 y'(t)\right)+g q_0-g z(t)+x'(t)^2+y'(t)^2+z'(t)^2\right)}{r^2-2 r x(t) \cos (\omega t)-2 a y(t) \sin (\omega t)+q_0^2-2 q_0 z(t)+x(t)^2+y(t)^2+z(t)^2}\\ y''(t)=\frac{(r \sin (\omega t)-y(t)) \left(r \omega \sin (\omega t) \left(\omega y(t)+2 x'(t)\right)+r \omega \cos (\omega t) \left(\omega x(t)-2 y'(t)\right)+g q_0-g z(t)+x'(t)^2+y'(t)^2+z'(t)^2\right)}{r^2-2 r x(t) \cos (\omega t)-2 r y(t) \sin (\omega t)+q_0^2-2 q_0 z(t)+x(t)^2+y(t)^2+z(t)^2}\\ z''(t)= \frac{-g r^2+r x(t) \cos (\omega t) \left(2 g+\omega ^2 q_0-\omega ^2 z(t)\right)+r y(t) \sin (\omega t) \left(2 g+\omega ^2 q_0-\omega ^2 z(t)\right)+(q_0-z(t)) \left(2 r \omega \sin (\omega t) x'(t)-2 r \omega \cos (\omega t) y'(t)+x'(t)^2+y'(t)^2+z'(t)^2\right)-g x(t)^2-g y(t)^2}{r^2-2 a x(t) \cos (\omega t)-2 r y(t) \sin (\omega t)+q_0^2-2 q_0 z(t)+x(t)^2+y(t)^2+z(t)^2} } $$
NOTA
La modelación se puede hacer considerando coordenadas polares como sugiere @mjqxxxx
Realizando
$$ p = q-l(\cos(\theta)\sin(\phi),\cos(\theta)\cos(\phi),\sin(\theta))\ \ \ \ \ (2) $$
y considerando ahora, el Lagrangiano
$$ L(p(\theta,\phi), \dot p(\theta,\phi),t) = \frac 12m \|\dot p\|^2-m g p\cdot\hat k $$
entonces las ecuaciones de movimiento resultantes son mucho más simples
$$ \cases{ \theta ''(t)= \frac{g \cos (\theta (t))-l \phi '(t)^2 \sin (\theta (t)) \cos (\theta (t))+\omega ^2 r \sin (\theta (t)) \sin (\omega t+\phi (t))}{l}\\ \phi ''(t)= 2 \phi '(t) \theta '(t) \tan (\theta (t))-\frac{\omega ^2 r \sec (\theta (t)) \cos (\omega t+\phi (t))}{l} } $$
Incluido un script de MATHEMATICA con animación
parms = {q0 -> 2, l -> 2, r -> 1, m -> 1, g -> 9.81, omega -> Pi/2};
ex = {1, 0, 0};
ey = {0, 1, 0};
ez = {0, 0, 1};
q = r (Cos[omega t] ex + Sin[omega t] ey) + q0 ez;
p = {x[t], y[t], z[t]};
L = 1/2 m D[p, t] . D[p, t] - m g p . ez + lambda ((p - q) . (p - q) - l^2);
mov = D[Grad[L, D[p, t]], t] - Grad[L, p];
equc = D[(p - q) . (p - q), t, t];
equs = Join[mov, {equc}];
sol = Solve[equs == 0, Join[D[p, t, t], {lambda}]][[1]];
equsmov = Thread[D[p, t, t] == (D[p, t, t] /. sol /. parms)];
cinits = {x[0] == r, y[0] == x'[0] == y'[0] == z'[0] == 0, z[0] == q0 - l} /. parms;
tmax = 10 omega /. parms;
solp = NDSolve[Join[equsmov, cinits], p, {t, 0, tmax}];
gr0 = ParametricPlot3D[Evaluate[p /. solp], {t, 0, tmax}, ViewPoint -> {4, 0, 2}, PlotStyle -> Thin]
Q[t_] := r (Cos[omega t] ex + Sin[omega t] ey) + q0 ez /. parms
X[t_] := Evaluate[x[t] /. solp]
Y[t_] := Evaluate[y[t] /. solp]
Z[t_] := Evaluate[z[t] /. solp]
d = 3;
dt = 0.05;
ListAnimate[Table[Show[
Graphics3D[{{Red, Sphere[{X[t][[1]], Y[t][[1]], Z[t][[1]]}, 0.1]},
{Line[{{0, 0, (q0 /. parms)}, {0, 0, 0}}]},
{Line[{{0, 0, (q0 /. parms)}, Q[t]}]},
{Line[{Q[t], {X[t][[1]], Y[t][[1]], Z[t][[1]]}}]}}],
PlotRange -> {{-d, d}, {-d, d}, {-1/2, d - 1/2}},
Boxed -> False,
Axes -> False,
ViewPoint -> {4, 0, 2}], {t, 0, tmax, dt}]
]