El modelo dinámico completo de un manipulador de $n$ GDL expresado en el espacio articular viene dado por:
$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{V}(\mathbf{q}, \dot{\mathbf{q}}) + \mathbf{G}(\mathbf{q}) = \boldsymbol{\tau}$$donde:
En la ecuación de la articulación $i$:
$$M_{ii}(\mathbf{q})\ddot{q}_i + \sum_{j \neq i} M_{ij}(\mathbf{q})\ddot{q}_j + V_i(\mathbf{q}, \dot{\mathbf{q}}) + G_i(\mathbf{q}) = \tau_i$$Se aísla la aceleración propia $\ddot{q}_i$:
$$\ddot{q}_i = \frac{1}{M_{ii}(\mathbf{q})} \left[ \tau_i - \tau_{\text{dist}, i} \right]$$donde el término de perturbación acoplada engloba la inercia cruzada, Coriolis y gravedad:
$$\tau_{\text{dist}, i} = \sum_{j \neq i} M_{ij}(\mathbf{q})\ddot{q}_j + V_i(\mathbf{q}, \dot{\mathbf{q}}) + G_i(\mathbf{q})$$La planta desacoplada para el eje $i$ tiene la función de transferencia:
$$G_i(s) = \frac{Q_i(s)}{\mathcal{T}_i(s)} = \frac{1}{J s^2}$$(O bien $G_i(s) = \frac{1}{s(J s + B)}$ si el enunciado incluye rozamiento viscoso $B$).
| Especificación dada | Fórmula de conversión | Comentario / Uso |
|---|---|---|
| $\omega_n$ y $\xi$ directos | Directo (ej: $\omega_n = 3\text{ rad/s}, \xi = 1.0$) | Caso más frecuente en exámenes recientes. |
| Sobreoscilación $M_p$ | $$\xi = \frac{-\ln(M_p)}{\sqrt{\pi^2 + \ln^2(M_p)}}$$ | Para $M_p = 5\% = 0.05 \implies \xi \approx 0.69$. Para $M_p = 0 \implies \xi \ge 1.0$. |
| Tiempo establecimiento $t_s$ (al 2%) | $$\omega_n = \frac{4}{\xi t_s}$$ | Frecuencia natural del par de polos complejos dominantes. |
| Tiempo de pico $t_p$ | $$\omega_n = \frac{\pi}{t_p \sqrt{1-\xi^2}}$$ | Instante en el que se alcanza el máximo transitorio. |
Como el PID añade un integrador $\frac{1}{s}$, la ecuación de lazo cerrado es de tercer orden. Se ubica el par dominante en $s^2 + 2\xi\omega_n s + \omega_n^2 = 0$ y un tercer polo real rápido en $s = -p_o$:
Multiplicando ambos factores:
$$P_{\text{des}}(s) = (s^2 + 2\xi\omega_n s + \omega_n^2)(s + p_o) = \mathbf{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}$$El controlador PID en su forma clásica interactúa con la planta:
$$C_{\text{PID}}(s) = K_p + K_d s + \frac{K_i}{s} = \frac{K_d s^2 + K_p s + K_i}{s}$$La ecuación característica del lazo cerrado es $1 + C_{\text{PID}}(s) G_i(s) = 0$:
$$s^3 + \frac{K_d + B}{J} s^2 + \frac{K_p}{J} s + \frac{K_i}{J} = 0$$Si el motor entrega par proporcional a la intensidad $\tau = K_t I$ (como en el examen de Junio 2018 P2):
Si un eje soporta par de gravedad dependiente de la posición ($G_i(q) \neq 0$), exigirle al término integral $K_i$ que levante todo el peso estático causa dos problemas:
La solución óptima de examen es añadir una compensación estática en lazo abierto:
$$\tau_i = \tau_{\text{PID}, i} + \hat{G}_i(q)$$donde $\hat{G}_i(q) = G_i(q)$ cancela de forma inmediata el par de gravedad, dejando que el lazo cerrado PID actúe exclusivamente sobre la inercia $\ddot{q}$.
Cuando el robot es de acoplamiento directo y se desea seguimiento riguroso de trayectorias a alta velocidad, se utiliza la **Linealización por Realimentación Exacta**:
Sustituyendo esta ley en la dinámica del manipulador:
$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}} + \mathbf{V} + \mathbf{G} = \mathbf{M}(\mathbf{q})\left[ \ddot{\mathbf{q}}_d + \mathbf{K}_v \dot{\mathbf{e}} + \mathbf{K}_p \mathbf{e} \right] + \mathbf{V} + \mathbf{G}$$Como $\mathbf{M}(\mathbf{q})$ es invertible $\forall \mathbf{q}$:
$$\mathbf{\ddot{e} + K_v \dot{e} + K_p e = 0}$$Para fijar polos con $\omega_{ni}$ y $\xi_i$ en cada articulación, las matrices de ganancia son estrictamente diagonales:
$$\mathbf{K}_p = \begin{bmatrix} \omega_{n1}^2 & 0 \\ 0 & \omega_{n2}^2 \end{bmatrix}, \quad \mathbf{K}_v = \begin{bmatrix} 2\xi_1\omega_{n1} & 0 \\ 0 & 2\xi_2\omega_{n2} \end{bmatrix}$$
Todo este procedimiento está programado en la función de MATLAB:
controlador_desde_modelo_dinamico.m
% EJEMPLO DE USO EN MATLAB (Octubre 2025):
syms q1 q2 real;
M = [1.0 + 3.0*cos(q2), 0.5 + 1.5*cos(q2);
0.5 + 1.5*cos(q2), 2.5];
G = [0; 2*9.81*sin(q2)];
% Actuadores directos con Kt = 10 N*m/A:
act = struct('Kt', 10.0, 'N', 1.0, 'Bm', 2e-5);
% Especificaciones: wn=3 rad/s, xi=1.0, po=4*wn:
specs = struct('wn', 3.0, 'xi', 1.0, 'po_factor', 4.0);
% Cálculo de ganancias:
resultado = controlador_desde_modelo_dinamico(M, G, act, specs);