ROBÓTICA — EXAMEN TERCERA CONVOCATORIA
Resolución Oficial Detallada: Problema 3 (Dinámica y Control PID)
Grado en Ingeniería de Tecnologías Industriales
Universidad de Sevilla | 29 de Octubre de 2025
Calificación: 3.5 Puntos (Total Examen: 10)
ENUNCIADO OFICIAL DEL EXAMEN
1.- (2.0 puntos). Dado el sistema descrito por las ecuaciones dinámicas siguientes: $$\begin{aligned} (m_1 + 2m_2\cos q_2)\ddot{q}_1 + (m_3 + m_2\cos q_2)\ddot{q}_2 - m_2\sin q_2\,\dot{q}_1\dot{q}_2 - m_2\sin q_2(\dot{q}_1 + \dot{q}_2)\dot{q}_2 &= T_1 \\ (m_3 + m_2\cos q_2)\ddot{q}_1 + m_4\ddot{q}_2 - m_2\sin q_2\,\dot{q}_1^2 + 2g\sin q_2 &= T_2 \end{aligned}$$ Donde $m_i$ son constantes, $g$ es la gravedad, y $q_i$, $T_i$ son las coordenadas generalizadas y las fuerzas/pares de articulación.
a) Escribir el modelo matricialmente: Determinar la matriz de inercia, el vector de fuerzas de Coriolis y centrífugas, y el vector de fuerzas gravitatorias.
b) Determinar los modelos aproximados lineales de cada articulación.
c) Diseñar un control PID para dichos modelos, para el caso de $m_1 = 1\text{ kg}$, $m_2 = 1.5\text{ kg}$, $m_3 = 2\text{ kg}$, $m_4 = 2.5\text{ kg}$.

2.- (1.5 puntos). Sea el robot de la figura adjunta. Calcular los parámetros dinámicos del eslabón 2 que se tendrían que introducir en el algoritmo de Newton-Euler utilizado en clase para poder obtener un modelo dinámico del robot. Para ello, considérese que los eslabones están construidos con varillas macizas de sección circular de radio $R = 0.05\text{ m}$ y densidad lineal $\rho = 1.5\text{ kg/m}$.

PROBLEMA 1: MODELADO MATRICIAL, LINEALIZACIÓN Y CONTROL PID

1.1. Formulación Matricial del Modelo Dinámico

La ecuación canónica del modelo dinámico de un robot manipulador en coordenadas articulares viene dada por: $$M(q)\ddot{q} + V(q, \dot{q}) + G(q) = T \quad \iff \quad M(q)\ddot{q} + C(q, \dot{q})\dot{q} + G(q) = T$$ donde $q = [q_1, q_2]^T$ y $T = [T_1, T_2]^T$.

1. Matriz de Inercia $M(q)$

Agrupando los coeficientes que multiplican a $\ddot{q}_1$ y $\ddot{q}_2$:

$$M(q) = \begin{bmatrix} m_1 + 2m_2\cos q_2 & m_3 + m_2\cos q_2 \\ m_3 + m_2\cos q_2 & m_4 \end{bmatrix}$$
✓ Propiedad fundamental: $M(q) = M^T(q)$ es simétrica y estrictamente definida positiva para configuraciones admisibles.
2. Términos de Velocidad $V(q, \dot{q})$ y Gravedad $G(q)$

Reescribiendo los términos de velocidad de la ecuación 1:

$$-m_2\sin q_2\dot{q}_1\dot{q}_2 - m_2\sin q_2(\dot{q}_1+\dot{q}_2)\dot{q}_2 = -m_2\sin q_2(2\dot{q}_1\dot{q}_2 + \dot{q}_2^2)$$ $$V(q, \dot{q}) = \begin{bmatrix} -m_2\sin q_2(2\dot{q}_1\dot{q}_2 + \dot{q}_2^2) \\ -m_2\sin q_2\,\dot{q}_1^2 \end{bmatrix}$$ $$G(q) = \begin{bmatrix} 0 \\ 2g\sin q_2 \end{bmatrix}$$
Forma matricial alternativa mediante matriz de Coriolis $C(q, \dot{q})$: $$C(q, \dot{q}) = \begin{bmatrix} -m_2\sin q_2\,\dot{q}_2 & -m_2\sin q_2(\dot{q}_1 + \dot{q}_2) \\ m_2\sin q_2\,\dot{q}_1 & 0 \end{bmatrix} \implies C(q,\dot{q})\dot{q} = V(q, \dot{q})$$ Con esta formulación, la matriz $\dot{M}(q) - 2C(q, \dot{q})$ es antisimétrica ($\dot{M}-2C = -(\dot{M}-2C)^T$), garantizando la pasividad del sistema.

1.2. Modelos Aproximados Lineales de cada Articulación (Desacoplo)

En el control articular industrial descentralizado (SISO), se desprecia el acoplamiento cruzado de aceleraciones ($\ddot{q}_j$ sobre el eje $i$) y las no linealidades de velocidad (Coriolis y centrífugas) al asumirse velocidades de régimen moderadas o compensadas por prealimentación:

Eje 1: Inercia Variable con $q_2$ (Sin Gravedad)

La inercia aparente vista por el motor 1 es:

$$J_1(q_2) = m_1 + 2m_2\cos q_2$$ Para un diseño robusto y conservador, se adopta la inercia máxima (caso peor), garantizando que el bucle no desestabilice: $$J_{1,\max} = m_1 + 2m_2 = 1.0 + 2(1.5) = 4.0\text{ kg}\cdot\text{m}^2$$ $$\frac{Q_1(s)}{T_1(s)} = \frac{1}{J_{1,\max} s^2} = \frac{1}{4.0\, s^2}$$
Eje 2: Inercia Constante + Gravedad

La inercia propia del eje 2 es constante:

$$J_2 = m_4 = 2.5\text{ kg}\cdot\text{m}^2$$ El término gravitatorio $G_2(q_2) = 2g\sin q_2$ se compensa activamente mediante prealimentación (*feedforward*): $T_2 = T_{2,\text{PID}} + 2g\sin q_2$. La dinámica vista por el PID es lineal: $$\frac{Q_2(s)}{T_{2,\text{PID}}(s)} = \frac{1}{J_2 s^2} = \frac{1}{2.5\, s^2}$$

1.3. Diseño y Sintonización del Controlador PID

Considerando un controlador PID estándar con acción filtrada/ideal $C(s) = K_p + K_d s + \frac{K_i}{s} = \frac{K_d s^2 + K_p s + K_i}{s}$: $$\text{Ecuación característica en lazo cerrado: } 1 + C(s)\frac{1}{J s^2} = 0 \iff s^3 + \frac{K_d}{J}s^2 + \frac{K_p}{J}s + \frac{K_i}{J} = 0$$ Se especifica un comportamiento en lazo cerrado con coeficiente de amortiguamiento crítico $\xi = 1$ (cero sobreoscilación), pulsación natural $\omega_n$ y un polo real rápido auxiliar en $p_o = 4\omega_n$: $$P_d(s) = (s^2 + 2\xi\omega_n s + \omega_n^2)(s + p_o) = s^3 + (2\xi\omega_n + p_o)s^2 + (\omega_n^2 + 2\xi\omega_n p_o)s + p_o\omega_n^2$$

Fórmulas de Sintonización Analítica Cerrada: $$K_d = J(2\xi\omega_n + p_o), \qquad K_p = J(\omega_n^2 + 2\xi\omega_n p_o), \qquad K_i = J(p_o\omega_n^2)$$
Articulación Inercia Nominal $J$ Diseño ($\xi, \omega_n, p_o$) $K_p$ [N·m/rad] $K_d$ [N·m·s/rad] $K_i$ [N·m/(rad·s)]
Eje 1 $J_{1,\max} = 4.0\text{ kg}\cdot\text{m}^2$ $\xi=1, \omega_{n1}=3\text{ rad/s}, p_{o1}=12$ $4.0 \cdot (9 + 72) = \mathbf{324.0}$ $4.0 \cdot (6 + 12) = \mathbf{72.0}$ $4.0 \cdot (12 \cdot 9) = \mathbf{432.0}$
Eje 2 $J_2 = 2.5\text{ kg}\cdot\text{m}^2$ $\xi=1, \omega_{n2}=4\text{ rad/s}, p_{o2}=16$ $2.5 \cdot (16 + 128) = \mathbf{360.0}$ $2.5 \cdot (8 + 16) = \mathbf{60.0}$ $2.5 \cdot (16 \cdot 16) = \mathbf{640.0}$

PROBLEMA 2: PARÁMETROS DINÁMICOS NEWTON-EULER DEL ESLABÓN 2 (MANIVELA)

Robot Manivela Octubre 2025 P3
Figura 1: Esquema cinemático y geométrico del robot con eslabón 2 tipo manivela (3 tramos cilíndricos de radio $R$ y densidad lineal $\rho$).

2.1. Descomposición Geométrica y Masa del Eslabón 2

El eslabón 2 está compuesto por 3 tramos de varilla maciza cilíndrica de sección circular de radio $R = 0.05\text{ m}$ y densidad lineal $\rho = 1.5\text{ kg/m}$:

Masa Total del Eslabón 2: $$m_2 = m_A + m_B + m_C = \rho\left(\frac{L_{2B}}{2}\right) + \rho L_{2A} + \rho\left(\frac{L_{2B}}{2}\right) = \mathbf{\rho(L_{2A} + L_{2B})}$$

2.2. Centro de Masas del Eslabón 2 (Vector $\mathbf{s}_{22}$)

Calculamos las coordenadas del centro de gravedad $C_2$ respecto al sistema de referencia del eslabón $\{2\}$:

Coordenada $X_{C2}$: $$X_{C2} = \frac{m_A\frac{L_{2B}}{4} + m_B\frac{L_{2B}}{2} + m_C\frac{3L_{2B}}{4}}{m_2} = \mathbf{\frac{L_{2B}}{2}}$$ (Resultado evidente por la perfecta simetría horizontal del eslabón manivela).
Coordenada $Y_{C2}$: $$Y_{C2} = \frac{m_A(0) + (\rho L_{2A})\frac{L_{2A}}{2} + \left(\rho\frac{L_{2B}}{2}\right)L_{2A}}{\rho(L_{2A} + L_{2B})} = \mathbf{\frac{L_{2A}\left(\frac{L_{2A}}{2} + \frac{L_{2B}}{2}\right)}{L_{2A} + L_{2B}} = \frac{L_{2A}}{2}\frac{L_{2A}+L_{2B}}{L_{2A}+L_{2B}} + \dots = \frac{L_{2A}(L_{2A} + L_{2B}/2)}{L_{2A} + L_{2B}}}$$

En el algoritmo de Newton-Euler (como el implementado en NE_R3GDL.m), el vector $\mathbf{s}_{22}$ representa la posición del centro de gravedad del eslabón 2 expresada en los ejes del sistema $\{2\}$: $$\mathbf{s}_{22} = \begin{bmatrix} X_{C2} - L_{2B} \\ Y_{C2} - L_{2A} \\ 0 \end{bmatrix} = \begin{bmatrix} -\frac{L_{2B}}{2} \\ -\frac{L_{2A}^2 / 2}{L_{2A} + L_{2B}} \\ 0 \end{bmatrix}$$

2.3. Tensor de Inercia Baricéntrico $\mathbf{I}_{22}$

Para un cilindro macizo de radio $R$, masa $m$ y longitud $L$: $$I_{\parallel} = \frac{1}{2}m R^2 \quad (\text{eje longitudinal}), \qquad I_{\perp} = \frac{1}{12}m L^2 + \frac{1}{4}m R^2 \quad (\text{ejes transversales})$$ Aplicando el Teorema de Steiner ($\mathbf{I}_{C} = \sum (\mathbf{I}_{i,C_i} + m_i [(\mathbf{d}_i^T\mathbf{d}_i)\mathbf{I}_3 - \mathbf{d}_i \mathbf{d}_i^T])$) respecto al centro de masas global $C_2$:

Tramo Masa $m_i$ $I_{xx, i}^{(0)}$ $I_{yy, i}^{(0)}$ $I_{zz, i}^{(0)}$ Vector distancia $\mathbf{d}_i = C_i - C_2$
Tramo A $\rho \frac{L_{2B}}{2}$ $\frac{1}{2}m_A R^2$ $\frac{1}{12}m_A (\frac{L_{2B}}{2})^2 + \frac{1}{4}m_A R^2$ $\frac{1}{12}m_A (\frac{L_{2B}}{2})^2 + \frac{1}{4}m_A R^2$ $[-\frac{L_{2B}}{4}, -Y_{C2}, 0]^T$
Tramo B $\rho L_{2A}$ $\frac{1}{12}m_B L_{2A}^2 + \frac{1}{4}m_B R^2$ $\frac{1}{2}m_B R^2$ $\frac{1}{12}m_B L_{2A}^2 + \frac{1}{4}m_B R^2$ $[0, \frac{L_{2A}}{2} - Y_{C2}, 0]^T$
Tramo C $\rho \frac{L_{2B}}{2}$ $\frac{1}{2}m_C R^2$ $\frac{1}{12}m_C (\frac{L_{2B}}{2})^2 + \frac{1}{4}m_C R^2$ $\frac{1}{12}m_C (\frac{L_{2B}}{2})^2 + \frac{1}{4}m_C R^2$ $[\frac{L_{2B}}{4}, L_{2A} - Y_{C2}, 0]^T$
Tensor de Inercia Final para Newton-Euler: $$I_{22} = \begin{bmatrix} I_{xx} & I_{xy} & 0 \\ I_{yx} & I_{yy} & 0 \\ 0 & 0 & I_{zz} \end{bmatrix}$$ Dado que los tramos A y C están dispuestos asimétricamente en el plano $XY$ respecto a las diagonales, el producto de inercia $I_{xy} = -\sum m_i d_{xi} d_{yi} \neq 0$, mientras que $I_{xz} = I_{yz} = 0$ debido a la planicidad exacta del eslabón en el plano $XY$.

2.4. Implementación y Verificación en MATLAB

%% Parámetros Newton-Euler calculados para el eslabón 2:
m2 = rho * (L2A + L2B);
s22 = [-L2B/2; yC - L2A; 0];
I22 = I_steiner_tramoA + I_steiner_tramoB + I_steiner_tramoC;
fprintf('Masa m2: %.3f kg\n', double(subs(m2, [rho, L2A, L2B], [1.5, 0.4, 0.3])));
% Listo para alimentar directamente el bucle cinemático y dinámico de NE_R3GDL.m