ROBÓTICA — PRIMERA CONVOCATORIA (JUNIO 2018)
Resolución Oficial: Problema 2 (Dinámica Robot Plano RP y Control PID)
Grado en Ingeniería de Tecnologías Industriales
Universidad de Sevilla | 5 de Junio de 2018
Tiempo: 1h 30m | Calificación: 4.0 Puntos
ENUNCIADO DEL PROBLEMA (4.0 PUNTOS)
Sea un robot plano de dos grados de libertad de accionamiento directo con articulaciones RP como el de la figura. Supóngase que las piezas tienen una densidad constante.
- Pieza A (Base cilíndrica): $^A M = 1\text{ kg}$, $^{A}I_{11} = {^{A}}I_{22} = 0.0308\text{ kg}\cdot\text{m}^2$, $^{A}I_{33} = 0.02\text{ kg}\cdot\text{m}^2$, $H_B = 0.25\text{ m}$, $H = 1\text{ m}$.
- Pieza B (Barra prismática): $^B M = 1\text{ kg}$, $^{B}I_{11} = {^{B}}I_{33} = 0.5\text{ kg}\cdot\text{m}^2$, $^{B}I_{22} = 0.00021\text{ kg}\cdot\text{m}^2$, longitud $2L = 1\text{ m}$ ($L = 0.5\text{ m}$).
- Motores de accionamiento directo: $J_{m1} = J_{m2} = 0.025\text{ kg}\cdot\text{m}^2$, $B_{m1} = B_{m2} = 2\cdot 10^{-5}\text{ Nm/(rad/s)}$, $K_{t1} = 10\text{ Nm/A}$, $K_{t2} = 10\text{ N/A}$.

Se pide:
1) (1.5 puntos). Calcular su modelo dinámico directo, tomando como señales de control las intensidades de los motores ($I_1, I_2$).
2) (1.0 punto). Programar la función ModeloDinamico_R.m y crear un simulador del robot en Simulink (robot_RP.mdl).
3) (1.5 puntos). Implementar un controlador tipo PID diseñado, y mostrar los resultados de simulación.
Robot RP Junio 2018
Figura 1: Vistas frontal y superior del manipulador plano R-P de accionamiento directo con indicación de parámetros.

1. MODELO DINÁMICO DIRECTO CON INTENSIDADES (1.5 PUNTOS)

1.1. Análisis Cinemático y Energético (Euler-Lagrange)

El robot se mueve en el plano horizontal $XY$ perpendicular al vector de gravedad terrestre $\mathbf{g} = [0, 0, -g]^T$. Por consiguiente, la energía potencial gravitatoria en las coordenadas articulares es nula: $$V(q_1, q_2) = 0 \implies G(q) = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

La posición del centro de masas de la barra deslizante B respecto al eje de giro base en función de $q_1$ (giro) y $q_2$ (traslación prismática) es: $$x_B = q_2 \cos q_1, \qquad y_B = q_2 \sin q_1$$ Derivando respecto al tiempo, la velocidad lineal del CoM de B es: $$v_B^2 = \dot{x}_B^2 + \dot{y}_B^2 = (\dot{q}_2 \cos q_1 - q_2 \sin q_1 \dot{q}_1)^2 + (\dot{q}_2 \sin q_1 + q_2 \cos q_1 \dot{q}_1)^2 = \dot{q}_2^2 + q_2^2 \dot{q}_1^2$$

1.2. Energía Cinética y Matriz de Inercia

La energía cinética total del sistema mecánico y los rotores de los motores es:

$$T = \frac{1}{2}\underbrace{(^{A}I_{33} + {^{B}}I_{33} + J_{m1})}_{J_0}\dot{q}_1^2 + \frac{1}{2} {^{B}}M (\dot{q}_2^2 + q_2^2 \dot{q}_1^2) + \frac{1}{2} J_{m2} \dot{q}_2^2$$ $$T = \frac{1}{2}\left( J_0 + {^{B}}M q_2^2 \right) \dot{q}_1^2 + \frac{1}{2}\left( {^{B}}M + J_{m2} \right) \dot{q}_2^2$$ donde $J_0 = 0.02 + 0.5 + 0.025 = 0.545\text{ kg}\cdot\text{m}^2$ y $M_{2,\text{eff}} = 1.0 + 0.025 = 1.025\text{ kg}$. $$\mathbf{M(q) = \begin{bmatrix} J_0 + {^{B}}M q_2^2 & 0 \\ 0 & {^{B}}M + J_{m2} \end{bmatrix} = \begin{bmatrix} 0.545 + q_2^2 & 0 \\ 0 & 1.025 \end{bmatrix}}$$

1.3. Términos de Coriolis y Centrífugos

$$\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{q}_1}\right) - \frac{\partial T}{\partial q_1} = (J_0 + {^{B}}M q_2^2)\ddot{q}_1 + 2 {^{B}}M q_2 \dot{q}_1 \dot{q}_2$$ $$\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{q}_2}\right) - \frac{\partial T}{\partial q_2} = ({^{B}}M + J_{m2})\ddot{q}_2 - {^{B}}M q_2 \dot{q}_1^2$$

1.4. Modelo en Espacio de Estados y Aceleración Directa

Incorporando el rozamiento viscoso de los motores y los pares generados por las corrientes de control ($\tau_1 = K_{t1} I_1 = 10 I_1$, $\tau_2 = K_{t2} I_2 = 10 I_2$):

$$\begin{aligned} (0.545 + q_2^2)\ddot{q}_1 + 2 q_2 \dot{q}_1 \dot{q}_2 + 2\cdot 10^{-5}\dot{q}_1 &= 10 I_1 \\ 1.025 \ddot{q}_2 - q_2 \dot{q}_1^2 + 2\cdot 10^{-5}\dot{q}_2 &= 10 I_2 \end{aligned}$$
Ecuaciones de Aceleración Directa para ModeloDinamico_R.m: $$\ddot{q}_1 = \frac{10 I_1 - 2 q_2 \dot{q}_1 \dot{q}_2 - 2\cdot 10^{-5}\dot{q}_1}{0.545 + q_2^2}$$ $$\ddot{q}_2 = \frac{10 I_2 + q_2 \dot{q}_1^2 - 2\cdot 10^{-5}\dot{q}_2}{1.025}$$

2. FUNCIÓN MATLAB ModeloDinamico_R.m Y SIMULINK (1.0 PUNTO)

function [qdd] = ModeloDinamico_R(q, qd, I)
  % MODELODINAMICO_R: Calcula las aceleraciones articulares qdd a partir del
  % estado [q, qd] y el vector de intensidades de motor I = [I1; I2].
  
  % Parámetros físicos del robot:
  J0     = 0.545;     % kg*m^2 (I_A33 + I_B33 + J_m1)
  M_B    = 1.0;       % kg
  M2_eff = 1.025;     % kg (M_B + J_m2)
  Bm1    = 2e-5;      % Nm/(rad/s)
  Bm2    = 2e-5;      % N/(m/s)
  Kt1    = 10.0;      % Nm/A
  Kt2    = 10.0;      % N/A
  
  q1  = q(1);  q2  = q(2);
  qd1 = qd(1); qd2 = qd(2);
  I1  = I(1);  I2  = I(2);
  
  % Dinámica articular directa:
  qdd1 = (Kt1 * I1 - 2 * M_B * q2 * qd1 * qd2 - Bm1 * qd1) / (J0 + M_B * q2^2);
  qdd2 = (Kt2 * I2 + M_B * q2 * (qd1^2) - Bm2 * qd2) / M2_eff;
  
  qdd = [qdd1; qdd2];
end

3. DISEÑO DEL CONTROLADOR PID Y RESULTADOS (1.5 PUNTOS)

3.1. Linealización Desacoplada y Peor Caso de Inercia

Eje 1 (Rotación): La inercia varía con la extensión $q_2 \in [0, L]$. Para garantizar estabilidad en todo el espacio de trabajo, se diseña para la máxima inercia (peor caso, $q_2 = 0.5\text{ m}$): $$J_{1,\max} = J_0 + {^{B}}M (0.5^2) = 0.545 + 0.25 = \mathbf{0.795\text{ kg}\cdot\text{m}^2}$$ $$\frac{Q_1(s)}{I_1(s)} = \frac{K_{t1}}{J_{1,\max} s^2} = \frac{10}{0.795 s^2} = \frac{12.578}{s^2}$$
Eje 2 (Prismático): La inercia es invariante: $$M_{2,\text{eff}} = 1.025\text{ kg}$$ $$\frac{Q_2(s)}{I_2(s)} = \frac{K_{t2}}{M_{2,\text{eff}} s^2} = \frac{10}{1.025 s^2} = \frac{9.756}{s^2}$$

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

Imponiendo amortiguamiento crítico $\xi = 1$ (cero sobreoscilación) y polo auxiliar rápido $p_o = 4\omega_n$: $$P_d(s) = (s + \omega_n)^2 (s + 4\omega_n) = s^3 + 6\omega_n s^2 + 9\omega_n^2 s + 4\omega_n^3$$ Comparando con el lazo cerrado de $G(s) = \frac{K_t}{J s^2}$ y PID $C(s) = \frac{K_d s^2 + K_p s + K_i}{s}$: $$K_d = \frac{J}{K_t} (6\omega_n), \qquad K_p = \frac{J}{K_t} (9\omega_n^2), \qquad K_i = \frac{J}{K_t} (4\omega_n^3)$$

Eje Inercia $J$ $\omega_n$ [rad/s] $K_p$ [A/rad o A/m] $K_d$ [A·s/rad o A·s/m] $K_i$ [A/(rad·s) o A/(m·s)]
Eje 1 ($q_1$) $0.795\text{ kg}\cdot\text{m}^2$ $\omega_{n1} = 4.0$ $\frac{0.795}{10} \cdot 9(16) = \mathbf{11.448}$ $\frac{0.795}{10} \cdot 6(4) = \mathbf{1.908}$ $\frac{0.795}{10} \cdot 4(64) = \mathbf{20.352}$
Eje 2 ($q_2$) $1.025\text{ kg}$ $\omega_{n2} = 5.0$ $\frac{1.025}{10} \cdot 9(25) = \mathbf{23.063}$ $\frac{1.025}{10} \cdot 6(5) = \mathbf{3.075}$ $\frac{1.025}{10} \cdot 4(125) = \mathbf{51.250}$

3.3. Gráficas de Simulación Temporal

Respuesta PID Robot RP
Figura 2: Evolución temporal de posición articular ($q_1, q_2$) ante escalones de referencia de 45° y 0.3 m en lazo cerrado no lineal.
Análisis de Resultados: El tiempo de establecimiento es inferior a $1.2\text{ segundos}$ en ambos ejes, con sobreoscilación prácticamente nula ($M_p < 2\%$) y error nulo en régimen permanente ($e_{ss} = 0$), demostrando la robustez del diseño frente a las fuerzas no lineales de Coriolis acopladas.