Solución Oficial y Detallada de Examen

Robótica — Grado en Ingeniería de Tecnologías Industriales (GITI) | Universidad de Sevilla
Convocatoria: 9 de Junio de 2025 | Problema B (3.5 Puntos)

Enunciado Oficial del Problema B

Sea un robot de dos grados de libertad de accionamiento directo con articulaciones RR como el de la figura. Supóngase que las piezas tienen una densidad constante.

Las piezas A y B son iguales en longitud (\(L_1=L_2=2L\)) con masa \({^AM} = 1.0\text{ Kg}\), tienen la misma densidad y momento de inercia que la B. Para dicho robot, y en función de los parámetros del mismo, se pide:

  1. (1.0 pts) Calcular su modelo dinámico directo, tomando como señales de control las intensidades de los motores. Se proporciona un fichero Calculo_NE.m como ayuda.
  2. (1.0 pts) Se proporciona un archivo de Simulink llamado robot_RR.mdl. Se pide programar la función ModeloDinamico_R.m y crear un simulador del robot propuesto utilizando las ecuaciones obtenidas del apartado anterior.
  3. (1.5 pts) Implementar un controlador tipo PID diseñado, y mostrar los resultados de simulación. Esto puede ser realizado mediante la ejecución del fichero graficas, siempre y cuando se respeten los nombres de las variables del espacio de trabajo de Matlab.
Esquema Técnico del Robot RR de Accionamiento Directo
Figura 1. Configuración cinemático-espacial del robot RR, detalle baricéntrico de la pieza B con ejes principales y tabla resumen de parámetros físicos y electromecánicos.

1. Deducción Rigurosa del Modelo Dinámico Directo (Apartado 1)

El manipulador es un robot espacial de 2 GDL con cinemática RR (Rotacional-Rotacional) de accionamiento directo (direct-drive, relación de transmisión \(r = 1\)):

Parámetros Numéricos Prescritos

Elemento Parámetro Físico Símbolo Valor Numérico Unidades
Eslabones A y B Longitud total (\(2L\)) \(L_1, L_2\) \(1.0\) \(\text{m}\)
Eslabones A y B Masa total \(m_1, m_2\) \(1.0\) \(\text{kg}\)
Eslabones A y B Distancia al centro de masas (\(L\)) \(l_{c1}, l_{c2}\) \(0.5\) \(\text{m}\)
Pieza B (y A) Inercia transversal baricéntrica \(I_{\perp} = {^B I_{11}} = {^B I_{33}}\) \(0.5\) \(\text{kg}\cdot\text{m}^2\)
Pieza B (y A) Inercia longitudinal baricéntrica \(I_{\parallel} = {^B I_{22}}\) \(0.00021\) \(\text{kg}\cdot\text{m}^2\)
Motores Direct-Drive Inercia de rotor \(J_{m1}, J_{m2}\) \(0.025\) \(\text{kg}\cdot\text{m}^2\)
Motores Direct-Drive Fricción viscosa de rotor \(B_{m1}, B_{m2}\) \(2.0 \times 10^{-5}\) \(\text{N}\cdot\text{m}/(\text{rad/s})\)
Motores Direct-Drive Constante de par electromagnético \(K_{t1}, K_{t2}\) \(10.0\) \(\text{N}\cdot\text{m/A}\)

Cinemática y Energía Cinética

1Eslabón 1 (Pieza A):

Al estar en el plano horizontal a altura \(z = 0\), su centro de masas se ubica en \(\mathbf{r}_{c1} = [l_{c1}\cos q_1, l_{c1}\sin q_1, 0]^T\). Su velocidad lineal es \(v_{c1}^2 = l_{c1}^2 \dot{q}_1^2\). Su velocidad angular es \(\boldsymbol{\omega}_1 = [0, 0, \dot{q}_1]^T\). El eje vertical es transversal a la barra, por lo que su momento de inercia respecto al eje de giro es \(I_{\perp} = 0.5\):

\[ K_1 = \frac{1}{2} m_1 l_{c1}^2 \dot{q}_1^2 + \frac{1}{2} I_{\perp} \dot{q}_1^2 + \frac{1}{2} J_{m1} \dot{q}_1^2 = \frac{1}{2} \left( m_1 l_{c1}^2 + I_{\perp} + J_{m1} \right) \dot{q}_1^2 = \frac{1}{2} (0.775) \dot{q}_1^2 \]

2Eslabón 2 (Pieza B):

La posición de su centro de masas en el espacio tridimensional es:

\[ \mathbf{r}_{c2} = \begin{bmatrix} (L_1 + l_{c2}\cos q_2)\cos q_1 \\ (L_1 + l_{c2}\cos q_2)\sin q_1 \\ l_{c2}\sin q_2 \end{bmatrix} \implies v_{c2}^2 = (L_1 + l_{c2}\cos q_2)^2 \dot{q}_1^2 + l_{c2}^2 \dot{q}_2^2 \]

La velocidad angular total de la barra 2 es \(\boldsymbol{\omega}_2 = \dot{q}_1 \hat{\mathbf{k}} + \dot{q}_2 \hat{\mathbf{e}}_2\). Proyectando sobre sus ejes baricéntricos principales:

Como \(I_{11} = I_{33} = I_{\perp}\), la energía cinética rotacional de la barra 2 es:

\[ K_{2,\text{rot}} = \frac{1}{2} I_{\parallel} (\dot{q}_1 \sin q_2)^2 + \frac{1}{2} I_{\perp} (\dot{q}_1^2 \cos^2 q_2 + \dot{q}_2^2) \]

Añadiendo la inercia del rotor del motor 2 montado en la articulación 2: \(\frac{1}{2} J_{m2} \dot{q}_2^2\).

Matriz de Inercia \(M(q)\)

Sumando todas las componentes, la energía cinética total no contiene términos cruzados \(\dot{q}_1 \dot{q}_2\), resultando una matriz de inercia diagonal:

\[ M(q) = \begin{bmatrix} M_{11}(q_2) & 0 \\ 0 & M_{22} \end{bmatrix} \]

Donde:

\[ \begin{aligned} M_{11}(q_2) &= m_1 l_{c1}^2 + I_{\perp} + J_{m1} + m_2 (L_1 + l_{c2}\cos q_2)^2 + I_{\perp}\cos^2 q_2 + I_{\parallel}\sin^2 q_2 \\ &= 0.775 + (1.0 + \cos q_2 + 0.25\cos^2 q_2) + 0.5\cos^2 q_2 + 0.00021\sin^2 q_2 \\ &= 1.775 + \cos q_2 + 0.75\cos^2 q_2 + 0.00021\sin^2 q_2 \quad [\text{kg}\cdot\text{m}^2] \\[6pt] M_{22} &= m_2 l_{c2}^2 + I_{\perp} + J_{m2} = 1.0(0.25) + 0.5 + 0.025 = 0.775 \quad [\text{kg}\cdot\text{m}^2] \end{aligned} \]

Fuerzas Centrífugas, de Coriolis y Gravitatorias

Derivando \(M_{11}(q_2)\) respecto a \(q_2\):

\[ \frac{\partial M_{11}}{\partial q_2} = -\sin q_2 - 1.5\cos q_2 \sin q_2 + 0.00042\sin q_2 \cos q_2 \approx -(1 + 1.49958\cos q_2)\sin q_2 \]

Aplicando las ecuaciones de Euler-Lagrange:

\[ c_1(q_2, \dot{q}_1, \dot{q}_2) = \frac{\partial M_{11}}{\partial q_2}\dot{q}_1 \dot{q}_2, \qquad c_2(q_2, \dot{q}_1) = -\frac{1}{2}\frac{\partial M_{11}}{\partial q_2}\dot{q}_1^2 \]

La gravedad actúa verticalmente (\(-\hat{\mathbf{k}}\)). El eslabón 1 se mueve horizontalmente (\(g_1 = 0\)). El eslabón 2 se eleva con altura de CoM \(z = l_{c2}\sin q_2\):

\[ U = m_2 g l_{c2}\sin q_2 \implies g_2(q_2) = \frac{\partial U}{\partial q_2} = m_2 g l_{c2}\cos q_2 = (1.0)(9.81)(0.5)\cos q_2 = 4.905\cos q_2 \quad [\text{N}\cdot\text{m}] \]

Ecuaciones en Función de las Corrientes \(i_1, i_2\)

En accionamiento directo (\(r=1\)), los pares entregados a las articulaciones son proporcionales a las corrientes eléctricas de control:

\[ \tau_1 = K_{t1} i_1 = 10 i_1, \qquad \tau_2 = K_{t2} i_2 = 10 i_2 \]

Incorporando la fricción viscosa \(B_{m1}\dot{q}_1\) y \(B_{m2}\dot{q}_2\), el modelo dinámico directo resultante es:

\[ \begin{aligned} \ddot{q}_1 &= \frac{1}{M_{11}(q_2)} \left[ 10\, i_1 - \left(\frac{\partial M_{11}}{\partial q_2}\right)\dot{q}_1 \dot{q}_2 - 2\cdot 10^{-5}\dot{q}_1 \right] \\[8pt] \ddot{q}_2 &= \frac{1}{0.775} \left[ 10\, i_2 + \frac{1}{2}\left(\frac{\partial M_{11}}{\partial q_2}\right)\dot{q}_1^2 - 4.905\cos q_2 - 2\cdot 10^{-5}\dot{q}_2 \right] \end{aligned} \]

2. Programación de la Función para Simulink (Apartado 2)

Para su inclusión en el bloque de MATLAB Function del archivo robot_RR.mdl, se ha creado el archivo oficial Matlab/03_Dinamica_y_Modelado/ModeloDinamico_R.m:

function [qdd] = ModeloDinamico_R(in) %% MODELO DINAMICO DIRECTO DEL ROBOT RR DE ACCIONAMIENTO DIRECTO % in = [ q1; q2; qd1; qd2; i1; i2 ] -> qdd = [ qdd1; qdd2 ] % 1. Desempaquetar variables de entrada q1 = in(1); q2 = in(2); qd1 = in(3); qd2 = in(4); i1 = in(5); i2 = in(6); % Intensidades de control (A) % 2. Parametros fisicos y constantes m2 = 1.0; lc2 = 0.5; g_acc = 9.81; I_perp = 0.5; I_long = 0.00021; Jm2 = 0.025; Bm1 = 2e-5; Bm2 = 2e-5; Kt1 = 10.0; Kt2 = 10.0; % 3. Pares aplicados por accionamiento directo (r = 1) tau1 = Kt1 * i1; tau2 = Kt2 * i2; % 4. Matriz de inercia M(q) M11 = 0.775 + (1.0 + cos(q2) + 0.25*cos(q2)^2) + 0.5*cos(q2)^2 + 0.00021*sin(q2)^2; M22 = m2 * lc2^2 + I_perp + Jm2; % = 0.775 M = [M11, 0; 0, M22]; % 5. Terminos Coriolis, Centrifugos y Friccion dM11_dq2 = -sin(q2) - 1.5*cos(q2)*sin(q2) + 0.00042*sin(q2)*cos(q2); c1 = dM11_dq2 * qd1 * qd2; c2 = -0.5 * dM11_dq2 * qd1^2; V = [c1; c2]; F_visc = [Bm1 * qd1; Bm2 * qd2]; % 6. Par gravitatorio (Eje 1 horizontal -> g1=0; Eje 2 vertical -> g2) g2 = m2 * g_acc * lc2 * cos(q2); G = [0; g2]; % 7. Aceleraciones articulares directas: qdd = M \ (Tau - V - F_visc - G) qdd = M \ ([tau1; tau2] - V - F_visc - G); end

3. Diseño del Controlador PID y Resultados de Simulación (Apartado 3)

Dado que la matriz de inercia es diagonal (\(M_{12} = M_{21} = 0\)), el sistema dinámico se encuentra naturalmente desacoplado a bajas velocidades. Para cada articulación \(k\), la planta linealizada se modela como un doble integrador accionado en corriente:

\[ J_{k,\text{nom}} \ddot{q}_k = K_{tk} i_k \implies \frac{q_k(s)}{i_k(s)} = \frac{K_{tk}}{J_{k,\text{nom}} s^2} \]

Se implementa una ley de control PID para cada eje en términos de la corriente del motor:

\[ i_k(t) = K_{pk} e_k(t) + K_{ik} \int_0^t e_k(\tau) d\tau - K_{dk} \dot{q}_k(t) \quad [+ \text{feedforward de gravedad para eje 2}] \]

Sintonización Analítica por Asignación de Polos

La ecuación característica de lazo cerrado deseada para un comportamiento críticamente amortiguado (\(\xi = 1.0\)) con frecuencia natural \(\omega_n\) y polo rápido auxiliar \(p_o = 4\omega_n\) es:

\[ s^3 + \frac{K_{tk} K_{dk}}{J_k} s^2 + \frac{K_{tk} K_{pk}}{J_k} s + \frac{K_{tk} K_{ik}}{J_k} = (s^2 + 2\xi\omega_n s + \omega_n^2)(s + p_o) = s^3 + 6\omega_n s^2 + 9\omega_n^2 s + 4\omega_n^3 \]

Igualando coeficientes, las ganancias de control en corriente son:

\[ K_{pk} = \frac{J_{k,\text{nom}}}{K_{tk}} (9\omega_n^2), \qquad K_{dk} = \frac{J_{k,\text{nom}}}{K_{tk}} (6\omega_n), \qquad K_{ik} = \frac{J_{k,\text{nom}}}{K_{tk}} (4\omega_n^3) \]

Ganancias Numéricas Obtenidas

Eje Articular Inercia Nominal \(J_{k,\text{nom}}\) Frecuencia \(\omega_{nk}\) \(K_{pk}\) (A/rad) \(K_{dk}\) (A\(\cdot\)s/rad) \(K_{ik}\) (A/(rad\(\cdot\)s)) Compensación Gravitatoria
Eje 1 (Giro Base) \(2.525\text{ kg}\cdot\text{m}^2\) \(3.0\text{ rad/s}\) \(20.4525\) \(4.5450\) \(27.2700\) Ninguna (\(g_1 = 0\))
Eje 2 (Elevación) \(0.775\text{ kg}\cdot\text{m}^2\) \(4.0\text{ rad/s}\) \(11.1600\) \(1.8600\) \(19.8400\) \(i_{2,\text{grav}} = \frac{4.905\cos q_2}{10.0} = 0.4905\cos q_2\) A

Curvas de Simulación Temporal

Simulando la respuesta ante escalones de referencia \([q_{1,\text{ref}}, q_{2,\text{ref}}] = [45^\circ, 30^\circ]\) partiendo del reposo:

Curvas de Simulación PID
Figura 2. Evolución temporal de las variables articulares \(q_1(t), q_2(t)\), errores de seguimiento \(e_1(t), e_2(t)\), y señales de control de intensidad \(i_1(t), i_2(t)\).
Análisis de los Resultados de Simulación:
Robótica (GITI) — Universidad de Sevilla | Escuela Técnica Superior de Ingeniería