Del avión a la matriz A

El mapa: qué estamos construyendo

Todo el proyecto se resume en una sola línea. Un modelo de simulación de vuelo es una función que, dado el estado actual del avión y lo que hacen los mandos, te dice cómo cambia el estado:

• \(\mathbf x\) es el estado: el mínimo de números que describe al avión en un instante (dónde está, hacia dónde mira, qué tan rápido va y gira).
• \(\mathbf u\) son las entradas: lo que el piloto o el controlador puede mandar (flaperones, deflector de la deriva, aceleradores de los 4 motores).
• \(f\) es la física: aerodinámica + propulsión + gravedad, pasadas por las leyes de Newton y Euler y por la cinemática.

Si sabes evaluar \(f\), puedes simular (integrar en el tiempo), buscar el equilibrio (trim, donde \(\dot{\mathbf x}=0\)) y linealizar (sacar las matrices \(A\) y \(B\) para diseñar el control). Por eso todo el cuaderno gira alrededor de cómo construir \(f\) bien.

El vector de estado de 13

Son 3 + 3 + 4 + 3 = 13. Con ángulos de Euler (\(\phi,\theta,\psi\)) en vez del cuaternión serían 12. Para un tail-sitter usamos el cuaternión porque en vuelo estacionario (hover) el morro apunta hacia arriba, \(\theta = 90^\circ\), y ahí los ángulos de Euler fallan (capítulo 2 y 5).

Las entradas de tu avión

\(\delta_f\): los dos flaperones juntos (los dos abajo o los dos arriba). \(\delta_a\): los flaperones en diferencial (alerón, para alabear). \(\delta_r\): el deflector de la deriva. \(t_1\dots t_4\): acelerador de cada motor (ala derecha, ala izquierda, punta superior de la deriva, punta inferior).

Los bloques y cómo se conectan

Fíjate en dos cosas. La propulsión no solo empuja: su estela cambia el aire que ven el ala y la deriva (por eso hay una flecha de propulsión a aerodinámica). Y Newton y Euler no dan posiciones: dan aceleraciones. La cinemática las convierte en cambios de posición y orientación, y la integración cierra el lazo.

Tu notación y la del curso

En las clases FD02 a FD05 (que siguen a Cook, Flight Dynamics Principles), las velocidades totales en ejes cuerpo se llaman \(U, V, W\) y las minúsculas \(u, v, w\) son perturbaciones alrededor del equilibrio. En este cuaderno, como en Stevens & Lewis y en Beard & McLain, \(u, v, w\) son las velocidades totales. Cuando linealicemos (capítulo 10) usaremos \(\Delta u\) para la perturbación. Fuerzas \(X, Y, Z\) y momentos \(L, M, N\) son iguales en ambos.

Qué archivo hace qué

Marcos de referencia y rotaciones

Todos los errores de signo de la mecánica de vuelo vienen de mezclar marcos. Antes de cualquier ecuación hay que tener claro en qué ejes está escrita cada cantidad y cómo se pasa de unos a otros.

Los tres marcos que usarás

Los ángulos aerodinámicos salen directamente de las velocidades en ejes cuerpo:

Girar un vector de un marco a otro

Un mismo vector (por ejemplo, la gravedad) tiene componentes distintas en cada marco. La matriz de cosenos directores (DCM) \(C_{B/E}\) convierte componentes en ejes Tierra a componentes en ejes cuerpo:

Con ángulos de Euler, la orientación se construye con tres giros sucesivos, en este orden: primero guiñada \(\psi\) alrededor de \(z\), luego cabeceo \(\theta\) alrededor del nuevo \(y\), y por último alabeo \(\phi\) alrededor del nuevo \(x\). Es la secuencia 3-2-1.

La DCM es ortogonal: su inversa es su traspuesta. Para pasar de cuerpo a Tierra usas \(C_{E/B}=C_{B/E}^T\). Esto lo necesitarás para la posición.

El problema de los ángulos de Euler: gimbal lock

Con \(\theta=90^\circ\), el primer giro (\(\psi\), alrededor de la vertical) y el último (\(\phi\), alrededor del morro, que ahora también apunta vertical) giran alrededor del mismo eje. Se pierde un grado de libertad: hay infinitas combinaciones de \(\phi\) y \(\psi\) para la misma orientación, y las ecuaciones de las derivadas de los ángulos se dividen entre \(\cos\theta=0\). Para un avión normal no importa porque nunca vuela a \(90^\circ\). Para un tail-sitter es su punto de operación en hover.

Cuaterniones: la solución

Cualquier orientación se puede alcanzar con un único giro de ángulo \(\eta\) alrededor de un eje unitario \(\mathbf e\) (teorema de Euler). El cuaternión guarda ese giro con cuatro números:

Cuatro números para tres grados de libertad: sobra uno, y por eso existe la restricción de norma unidad. A cambio, no hay ninguna orientación singular. La DCM en función del cuaternión es:

Y para pasar de ángulos de Euler (3-2-1) a cuaternión, con \(c_\phi=\cos(\phi/2)\), \(s_\phi=\sin(\phi/2)\), etc.:

Newton: ecuaciones de traslación

Aquí sale la ecuación que preguntabas: cómo cambia la velocidad \(u\) (y \(v\), \(w\)). La idea entera cabe en una frase: Newton vale en ejes Tierra, pero queremos escribir la velocidad en ejes cuerpo, que giran. Todo lo demás es hacer la cuenta con cuidado.

Hipótesis

• Cuerpo rígido, masa constante (sin consumo de combustible: tu batería no cambia de masa).
• Tierra plana y sin rotación: el marco NED es inercial. Vale para vuelos cortos y lentos como los tuyos.
• Origen de los ejes cuerpo en el CG.

Paso 1: la derivada de un vector que gira

Escribimos la velocidad con sus componentes en ejes cuerpo: \(\mathbf V=u\,\hat{\mathbf i}+v\,\hat{\mathbf j}+w\,\hat{\mathbf k}\). Pero \(\hat{\mathbf i},\hat{\mathbf j},\hat{\mathbf k}\) son los ejes del avión, y giran con él a velocidad angular \(\boldsymbol\omega=(p,q,r)\).

Paso 2: regla del producto

Paso 3: desarrollar el producto vectorial

Paso 4: las tres ecuaciones

La ecuación de \(\dot u\), término a término

• \(rv-qw\): no es una fuerza. Aparece porque los ejes donde medimos \(u\) están girando. Si el avión cabecea hacia arriba (\(q>0\)) con \(w>0\), el eje \(x\) se aleja del vector velocidad: \(u\) disminuye y \(w\) aumenta aunque nada empuje. El término \(-qw\) lo recoge.
• \(X_{\text{aero}}\): componente \(x\) cuerpo de sustentación y resistencia: \(L\sin\alpha-D\cos\alpha\) por tira, como en tu wingStripTheory.m.
• \(\sum T_i\): los 4 empujes, todos por \(+x\).
• \(-g\sin\theta\): con el morro arriba, la gravedad frena. En hover (\(\theta=90^\circ\)) es \(-g\) entero.

Euler: ecuaciones de rotación

Lo mismo que en el capítulo 3, pero con giros. La versión de rotación de Newton dice que el momento total es la derivada del momento angular, otra vez en ejes Tierra. Y otra vez escribimos todo en ejes cuerpo y usamos el teorema de transporte.

Momento angular y tensor de inercia

Cada partícula \(dm\) a una distancia \(\mathbf r\) del CG se mueve con \(\boldsymbol\omega\times\mathbf r\). Sumando \(\mathbf r\times(\boldsymbol\omega\times\mathbf r)\,dm\) sobre todo el avión:

En ejes cuerpo, \(I\) es constante (el avión no cambia de forma). Esa es la razón principal para trabajar en ejes cuerpo. Si el avión es simétrico respecto al plano \(x\)-\(z\), entonces \(I_{xy}=I_{yz}=0\) y solo queda \(I_{xz}\).

Paso 1: teorema de transporte para \(\mathbf H\)

Paso 2: desarrollar las componentes

Paso 3: despejar \(\dot p\) y \(\dot r\)

\(\dot p\) y \(\dot r\) aparecen juntos en \(L\) y \(N\) por culpa de \(I_{xz}\). Es un sistema lineal de 2×2. La solución se escribe con los coeficientes \(\Gamma\) (Stevens & Lewis; Beard & McLain):

En código no hace falta escribirlos: wdot = I \ (M - cross(w, I*w)) resuelve el sistema directamente. Pero sí conviene que los sepas leer, porque así aparecen en los libros y en los exámenes.

Los rotores también tienen momento angular

Cada hélice con su rotor gira a \(\Omega_i\) (miles de rpm) alrededor de \(+x\). Su momento angular \(\mathbf h_i=s_iI_r\Omega_i\,\hat{\mathbf i}\) (con \(s_i=\pm1\) el sentido de giro) se suma al del avión:

Con \(\mathbf h=(H,0,0)\): \(-\boldsymbol\omega\times\mathbf h=(0,-rH,qH)\). Si cabeceas mientras las hélices giran, aparece un momento de guiñada. Si tus hélices giran en pares opuestos (como en tu defineHardware.m, spin = [1 1 -1 -1]), \(H\approx0\) y el efecto se cancela. Esto ya está en motorForces.m.

Cinemática y el modelo de 13 estados

Newton y Euler te dan \(\dot u,\dot v,\dot w,\dot p,\dot q,\dot r\). Faltan dos piezas: cómo cambia la orientación con \(p, q, r\) y cómo cambia la posición con \(u, v, w\). No hay física nueva aquí, solo geometría.

Velocidades de los ángulos de Euler

Ojo: \(p, q, r\) no son \(\dot\phi,\dot\theta,\dot\psi\). Las velocidades angulares están en ejes cuerpo; los ángulos de Euler se miden en ejes intermedios distintos.

Cinemática del cuaternión

En un intervalo pequeño \(dt\), el avión gira un ángulo \(|\boldsymbol\omega|dt\) alrededor de \(\boldsymbol\omega\). Componiendo ese giro pequeño con la orientación actual se llega a una ecuación lineal en \(\mathbf q\), sin senos ni divisiones:

Posición

La velocidad está en ejes cuerpo; la posición en ejes Tierra. Basta girar el vector con la traspuesta de la DCM:

El modelo completo

Fuerzas aerodinámicas

Las ecuaciones de los capítulos 3 a 5 valen para cualquier sólido: un avión, un dron o un ladrillo. Lo que hace que sea tu avión son las fuerzas y momentos. La aerodinámica tiene que entregar seis números en ejes cuerpo, en el CG:

Hay varias formas de construir \(\mathcal A\), de menos a más precisas: strip theory (lo que tienes), tablas de un método de paneles (VSPAERO, XFLR5), CFD o túnel de viento, y finalmente datos de vuelo. La función es la misma; cambia de dónde salen los números.

Nivel 1: strip theory, tira a tira

1. Velocidad local de cada tira en su cuarto de cuerda \(\mathbf r\): la del avión, más la debida al giro, más la de la estela de las hélices: \(\mathbf V_\ell=\mathbf V_B+\boldsymbol\omega\times\mathbf r+\Delta\mathbf V_{\text{estela}}\).
2. Ángulo de ataque local \(\alpha_\ell=\operatorname{atan2}(w_\ell,u_\ell)\) y presión dinámica \(\bar q_\ell=\tfrac12\rho|\mathbf V_\ell|^2\). La rotación \(p\) hace que el ala que baja vea más ángulo de ataque: eso es el amortiguamiento en alabeo.
3. Coeficientes de la sección: \(c_l=c_{l0}+a\,(\alpha_\ell+\tau\delta)\), con límite suave en \(c_{l,\max}\) para la entrada en pérdida, y \(c_d=c_{d0}+c_l^2/(\pi eA\!R)\).
4. Fuerzas de la tira: \(dL=\bar q_\ell c\,c_l\,dy\), \(dD=\bar q_\ell c\,c_d\,dy\), giradas a ejes cuerpo con \(\alpha_\ell\).
5. Momentos: \(\sum\mathbf r\times d\mathbf F\), más el momento propio del flap.

La pendiente de sustentación de un ala finita (Helmbold, incompresible) y la eficacia del flap (perfil delgado) son:

El flap hacia abajo da más sustentación, pero un momento de picado: la carga extra cae detrás del cuarto de cuerda.

Forma con coeficientes y derivadas de estabilidad

Los libros y las tablas de VSPAERO usan coeficientes adimensionales. Para linealizar o comparar con otros aviones se escriben así (lateral, por ejemplo):

Cada derivada mide cuánto cambia un coeficiente al mover una variable. Tus demos de MATLAB las calculan con diferencias centrales, igual que harás con la matriz \(A\) en el capítulo 10.

Nivel 2: tablas de VSPAERO o XFLR5

1. Modela el avión en OpenVSP: ala, deriva, fuselaje, con los flaperones como control surface sub-surfaces y grupos de control (flap: ganancias +1, +1; alerón: +1, −1).
2. Barre una rejilla: \(\alpha\) de −10° a 20°, \(\beta\) de −10° a 10°, y cada mando en 3 a 5 valores. Guarda \(C_X, C_Y, C_Z, C_l, C_m, C_n\) (o \(C_L, C_D, \dots\)) respecto a un punto de referencia.
3. Saca las derivadas de giro (\(C_{lp}\), \(C_{mq}\), \(C_{nr}\)) del análisis de estabilidad.
4. En MATLAB, carga la rejilla y usa griddedInterpolant para interpolar en vuelo.
5. Suma la resistencia viscosa (VSPAERO es no viscoso): cd0 de XFOIL o de la herramienta Parasite Drag de OpenVSP.

No hace falta recalcular para cada posición del CG. Si las tablas están referidas a un punto \(\mathbf r_{\text{ref}}\), el momento respecto al CG es:

Las derivadas de giro sí cambian con el CG, porque la velocidad de cada punto depende de dónde está el eje de giro.

Niveles 3 y 4

CFD (STAR-CCM+, OpenFOAM, SU2): úsalo donde VSPAERO falla, no en toda la envolvente. Por ejemplo, unos pocos casos a \(\alpha\) alto con la hélice girando (modelo de disco actuador o malla rotante). Túnel de viento: mide con hélices funcionando. Vuelo: identificación de sistemas (capítulo 13) corrige las derivadas que más importan para el control.

Propulsión: motor, hélice, estela

Este es el capítulo que más rápido fue en el chat, así que va despacio. La cadena física es:

El acelerador no manda empuje: manda voltaje. El motor acelera hasta que el par que da iguala al par que le pide la hélice. Esa velocidad de giro fija el empuje.

7.1 El motor sin escobillas (BLDC) desde cero

Un motor eléctrico son imanes y bobinas. Dos efectos lo explican todo:

• Hace de generador: al girar, las bobinas producen un voltaje proporcional a la velocidad, la fuerza contraelectromotriz \(e=K_e\,\Omega\).
• Hace de motor: la corriente por las bobinas produce un par proporcional a la corriente, \(Q_m=K_t\,(I-I_0)\). \(I_0\) es la corriente que se pierde en rozamientos y pérdidas magnéticas sin carga.

El circuito, en régimen permanente (se desprecia la inductancia), es una resistencia y un "generador":

El variador (ESC) hace PWM: con acelerador \(t\in[0,1]\), el motor ve en media \(V_m=t\,V_{\text{bat}}\). La batería cae de voltaje con la corriente total: \(V_{\text{bat}}=V_{oc}-R_b\sum I_i\). Por eso el empuje baja durante el vuelo.

7.2 La hélice: coeficientes adimensionales

Por análisis dimensional, empuje, par y potencia de una hélice de diámetro \(D\) que gira a \(n\) rev/s solo dependen de la relación de avance \(J\):

Las curvas \(C_T(J)\) y \(C_P(J)\) son lo que publica la base de datos de hélices de UIUC y los ficheros de APC. Tu loadPropData.m las lee. Si no las tienes, propBEMT.m las calcula a partir de la geometría de la pala.

7.3 Teoría de cantidad de movimiento (disco actuador)

El modelo más sencillo: la hélice es un disco de área \(A\) que acelera el aire. Lejos delante el aire va a \(V\), en el disco a \(V+v_i\), y lejos detrás a \(V+v_w\).

7.4 Teoría del elemento de pala (BEMT)

El disco actuador no sabe nada de la forma de la pala. BEMT corta la pala en anillos y trata cada trozo como un perfil 2D, igual que strip theory con el ala. Cada trozo a radio \(r\) ve una velocidad tangencial \(\Omega r\) y una axial \(V+v_i\). El ángulo que forman, \(\phi\), y el paso de la pala, \(\theta\), dan el ángulo de ataque de la sección \(\alpha=\theta-\phi\).

La \(v_i\) de cada anillo se busca para que el empuje del elemento de pala coincida con el de cantidad de movimiento del anillo (\(dT=4\pi r\rho F(V+v_i)v_i\,dr\), con \(F\) la pérdida de punta de Prandtl). Es justo lo que hace propBEMT.m.

7.5 Dónde se juntan motor y hélice: el punto de operación

El motor da un par que baja con la velocidad (porque la fuerza contraelectromotriz se come el voltaje). La hélice pide un par que sube con la velocidad (\(\propto\Omega^2\)). El rotor acelera o frena hasta que se cruzan:

7.6 Dinámica de las rpm: el motor no responde al instante

Fuera del equilibrio, la diferencia de pares acelera el rotor (inercia \(I_r\) de rotor + hélice):

Linealizando alrededor del punto de operación sale un sistema de primer orden con constante de tiempo:

Para motores de este tamaño, \(\tau_m\) suele estar entre 20 y 200 ms. Parece poco, pero en hover todo el control de actitud pasa por las rpm, así que \(\tau_m\) limita lo rápido que puede ser tu controlador (capítulo 11). Si metes \(\Omega_i\) como estados, el modelo pasa a 17 estados.

7.7 La estela (slipstream)

Detrás de la hélice el aire va más rápido y el chorro se estrecha. Con la teoría de cantidad de movimiento:

El chorro además gira en el sentido de la hélice (torbellino, swirl), porque el aire se lleva el par que le ha dado la hélice. Eso cambia el ángulo de ataque a un lado y otro del eje: un lado del ala gana sustentación y el otro la pierde.

7.8 Lo que los motores le hacen al avión

7.9 Usar los datos del fabricante

1. Motor (hoja de datos): \(K_V\), resistencia interna, \(I_0\) (con el voltaje de la prueba), corriente máxima.
2. ESC y batería: suma sus resistencias a \(R\). Usa el voltaje bajo carga.
3. Hélice: curvas \(C_T(J), C_P(J)\) de UIUC o APC, o BEMT con la geometría.
4. Tabla de banco del fabricante (acelerador, empuje, corriente, rpm): compara con checkMotorAgainstTest.m y ajusta kT.
5. Tu propio banco de empuje (capítulo 13): mide empuje, par, rpm y corriente. Es lo más fiable.

Ensamblar y simular \(f(\mathbf x,\mathbf u)\)

Ya tienes todas las piezas. Este es el contenido de tu Eqs_of_motion.m, todavía vacío. Escríbelo tú, mirando el capítulo 5; aquí tienes una versión para comparar después.

function xdot = tailsitter6dof(t, x, U, W, Fn, Mo, P)
% x = [xN yE zD  u v w  q0 q1 q2 q3  p q r]'   (13 estados)
V  = x(4:6);   qt = x(7:10);   w = x(11:13);
qt = qt / norm(qt);
C  = quat2dcm_be(qt);                 % ejes Tierra -> ejes cuerpo

% fuerzas y momentos (aero + motores + estela), en el CG, ejes cuerpo
F.u = V(1); F.v = V(2); F.w = V(3);
F.p = w(1); F.q = w(2); F.r = w(3);
F.rho = 1.225;                        % o atmósfera ISA con h = -x(3)
A  = tailsitterAero(W, Fn, Mo, F, U);
Fb = [A.Fx; A.Fy; A.Fz];
Mb = [A.Roll; A.Pitch; A.Yaw];        % ya incluye par de reacción y giroscópico

% Newton, Euler y cinemática (capítulo 5)
Vdot = Fb/P.m + C*[0; 0; P.g] - cross(w, V);
wdot = P.I \ (Mb - cross(w, P.I*w));
Om   = [0 -w(1) -w(2) -w(3); w(1) 0 w(3) -w(2); w(2) -w(3) 0 w(1); w(3) w(2) -w(1) 0];
qdot = 0.5*Om*qt + 0.5*(1 - qt'*qt)*qt;   % término que corrige la norma
pdot = C' * V;
xdot = [pdot; Vdot; qdot; wdot];
end

function C = quat2dcm_be(q)
q0=q(1); q1=q(2); q2=q(3); q3=q(4);
C = [q0^2+q1^2-q2^2-q3^2, 2*(q1*q2+q0*q3),     2*(q1*q3-q0*q2);
     2*(q1*q2-q0*q3),     q0^2-q1^2+q2^2-q3^2, 2*(q2*q3+q0*q1);
     2*(q1*q3+q0*q2),     2*(q2*q3-q0*q1),     q0^2-q1^2-q2^2+q3^2];
end

import numpy as np
from scipy.integrate import solve_ivp

def f(t, x, u, P, aero):
    V, q, w = x[3:6], x[6:10], x[10:13]
    q = q / np.linalg.norm(q)
    C = quat2dcm_be(q)
    Fb, Mb = aero(V, w, u)                  # tu modelo de fuerzas
    Vdot = Fb/P.m + C @ np.array([0, 0, P.g]) - np.cross(w, V)
    wdot = np.linalg.solve(P.I, Mb - np.cross(w, P.I @ w))
    p, qq, r = w
    Om = np.array([[0,-p,-qq,-r],[p,0,r,-qq],[qq,-r,0,p],[r,qq,-p,0]])
    qdot = 0.5*Om @ q + 0.5*(1 - q @ q)*q
    return np.concatenate([C.T @ V, Vdot, qdot, wdot])

sol = solve_ivp(f, (0, 10), x0, args=(u, P, aero), max_step=0.005)

Integrar en el tiempo

• ode45 (paso variable) es cómodo para explorar. Con dinámica de motores rápida el sistema se vuelve rígido; si va lento, prueba ode15s.
• Para probar controladores usa paso fijo (Runge-Kutta 4) al mismo ritmo que el controlador real o a un múltiplo. Así el simulador se parece al hardware: el controlador se ejecuta cada \(\Delta t\), con su retraso.
• Regla práctica: \(\Delta t\) al menos 10 veces menor que la constante de tiempo más rápida (la de los motores, \(\tau_m\approx\) 20–200 ms).

Pruebas de cordura (hazlas siempre)

Trim: buscar el equilibrio

Un punto de trim es un estado y unas entradas en los que el avión vuela en régimen permanente: velocidades y giros constantes. Matemáticamente, \(\dot{\mathbf x}=0\) para todos los estados menos la posición (que avanza) y, en un viraje, el rumbo.

Hover

Incógnitas: los 4 aceleradores, los mandos y la actitud. Con el avión simétrico, \(\theta=90^\circ\), los dos motores del ala iguales y los de la deriva iguales. Condiciones: \(\sum T_i=mg\) y momentos nulos. Los motores de arriba y abajo tienen que compensar cualquier momento de cabeceo de la estela sobre los flaps.

Crucero

Dada la velocidad \(V\) y vuelo nivelado: incógnitas \(\alpha\) (y por tanto \(\theta=\alpha\)), flaps \(\delta_f\) y acelerador. Condiciones: \(\dot u=\dot w=\dot q=0\). Tres ecuaciones, tres incógnitas.

Transición

Repite el trim para muchas velocidades entre 0 y la de crucero. La curva \(\theta(V)\) y los mandos de trim forman el "pasillo de transición". El controlador se programa a lo largo de esa curva (capítulo 11).

% incógnitas z = [alpha; df; thr]   (crucero nivelado a velocidad V)
res = @(z) trimResidual(z, V, W, Fn, Mo, P);   % devuelve [udot; wdot; qdot]
z0  = [deg2rad(4); 0; 0.5];
z   = fsolve(res, z0, optimoptions('fsolve','Display','iter'));

function r = trimResidual(z, V, W, Fn, Mo, P)
alpha = z(1);  U.df = z(2);  U.da = 0;  U.dr = 0;  U.thr = z(3)*ones(1,4);
x = stateFromTrim(V, alpha);          % u = V cos(alpha), w = V sin(alpha), theta = alpha
xd = tailsitter6dof(0, x, U, W, Fn, Mo, P);
r = xd([4 6 12]);                     % udot, wdot, qdot
end

Linealizar: las matrices \(A\) y \(B\)

Aquí llegas a la matriz que preguntabas. Cerca del trim, las pequeñas perturbaciones se comportan de forma lineal. Es una serie de Taylor de primer orden:

Cada casilla responde a una pregunta: si muevo un poco el estado \(j\), ¿cuánto cambia la derivada del estado \(i\)? Por ejemplo, \(\partial\dot u/\partial\theta=-g\cos\theta_0\): si levantas el morro, la gravedad empieza a frenar.

Cómo se calcula: diferencias centrales

function [A, B] = linearize(f, x0, u0)
% f(x,u) devuelve xdot. Diferencias centrales columna a columna.
n = numel(x0);  m = numel(u0);
A = zeros(n);   B = zeros(n, m);
for j = 1:n
    h = 1e-6*max(1, abs(x0(j)));  dx = zeros(n,1);  dx(j) = h;
    A(:,j) = (f(x0+dx, u0) - f(x0-dx, u0)) / (2*h);
end
for j = 1:m
    h = 1e-6*max(1, abs(u0(j)));  du = zeros(m,1);  du(j) = h;
    B(:,j) = (f(x0, u0+du) - f(x0, u0-du)) / (2*h);
end
end

Otra opción es derivar \(f\) simbólicamente (Symbolic Math Toolbox, jacobian; o SymPy en Python). Solo funciona si la aerodinámica es analítica. Con tablas o strip theory, numérico.

El problema del cuaternión: ¿13×13 o 12×12?

Si linealizas el modelo de 13 estados tal cual, \(A\) sale de 13×13 pero con una dirección "falsa": mover las cuatro componentes de \(\mathbf q\) a la vez cambia su norma, y eso no corresponde a ningún movimiento real. Aparece un autovalor espurio y el sistema no es controlable en esa dirección. Tres soluciones:

1. Crucero: linealiza con ángulos de Euler (12 estados). Lejos de \(\theta=90^\circ\) no hay problema. Es lo que hace el curso.
2. Error de actitud: escribe la actitud como trim por un giro pequeño, \(\mathbf q=\mathbf q_0\otimes[1,\ \tfrac12\delta\boldsymbol\theta]\), y usa los 3 componentes de \(\delta\boldsymbol\theta\) como estados. Queda 12×12 y vale en hover. Es lo que usan los filtros de Kalman de los autopilotos.
3. Marco de hover: define los ángulos de Euler respecto a un marco girado 90°. En hover esos ángulos son pequeños.

Qué haces con \(A\) y \(B\)

• eig(A) o damp(ss(A,B,C,D)): los modos. En crucero: corto periodo, fugoide, alabeo puro, balanceo holandés, espiral (clases FD13). En hover salen modos distintos, casi neutros y acoplados con la posición.
• rank(ctrb(A,B)): ¿puedes controlar todos los estados con tus mandos?
• lqr, place: diseño de controladores (clases 9–12 del curso).

Control para hardware real

Un controlador que funciona en simulación y vibra en el avión casi siempre falla por algo que el modelo no tenía: retrasos, respuesta lenta de los motores, saturación, ruido de los sensores. Este capítulo trata de lo que hace que el control funcione de verdad.

Arquitectura en cascada

Cada lazo es más rápido que el de fuera (típicamente de 3 a 5 veces). El más interno, el de velocidad angular, está limitado por la rapidez de los actuadores y por los retrasos. Así funcionan PX4 y ArduPilot.

Asignación de mandos (mixer)

Los lazos piden un empuje total y tres momentos. Hay que repartirlos entre tus 7 mandos. Linealizando la propulsión y la aerodinámica de los mandos:

\(B^+\) es la pseudoinversa (o una versión con pesos). La tabla de "eficacia de motores" que imprimía demo_tailsitter.m es justamente una parte de \(B_{\text{asig}}\). Ojo: depende de la velocidad. En hover los flaps solo funcionan por la estela; en crucero dominan sobre el empuje diferencial.

Métodos

• PID en cascada: simple, se ajusta en vuelo, es lo que usan los autopilotos abiertos.
• LQR sobre el modelo lineal: K = lqr(A,B,Q,R). Da un buen punto de partida para las ganancias. Clase 12 del curso.
• Programación de ganancias (gain scheduling): un controlador por cada punto del pasillo de transición, interpolando con la velocidad.
• INDI (inversión dinámica no lineal incremental): usa la aceleración angular medida y la eficacia de los mandos para calcular el incremento de mando. Depende poco del modelo aerodinámico y se ha usado mucho en investigación con tail-sitters. Es un buen tema para leer cuando domines lo básico.

Lista para que funcione en hardware

• Saturaciones de mandos y anti-windup en los integradores.
• Controlador discreto al ritmo real (250–1000 Hz en el lazo interno), con el retraso medido.
• Márgenes de estabilidad con modelo de actuador y retraso: como referencia típica, al menos 6 dB de ganancia y 45° de fase.
• Montecarlo: repite las simulaciones con derivadas, inercias, \(\tau_m\) y retrasos variando ±20–30 %. El controlador tiene que funcionar en todas.
• Filtros: pasa-bajo en la derivada y notch para la vibración de las hélices (capítulo 12).

Sensores, filtros y estimación

El controlador no ve el estado \(\mathbf x\): ve lo que miden los sensores, con ruido, sesgo y retraso. Un modelo para diseñar control que funcione en hardware tiene que incluir los sensores y el estimador.

Modelo de sensor

Medida = valor real + sesgo que deriva lentamente + ruido blanco. Añade la frecuencia de muestreo, la cuantización y el retraso. Los valores salen de la hoja de datos del sensor (densidad de ruido, inestabilidad del sesgo) o de una grabación con el avión quieto.

Vibración y filtros

Las hélices excitan la estructura a la frecuencia de paso de pala, \(f=B\cdot\text{rpm}/60\). Con 2 palas a 7000 rpm son 233 Hz. Si el giróscopo muestrea mal, esa vibración aparece como ruido de baja frecuencia (aliasing) y entra al lazo de control. Soluciones: aislar la IMU, filtro anti-aliasing, y un filtro notch que sigue las rpm (los autopilotos lo llaman harmonic notch).

Filtro complementario: la idea de toda la estimación

El giróscopo es bueno a corto plazo (integrado da el ángulo, pero deriva). El acelerómetro es bueno a largo plazo (ve la gravedad, pero es ruidoso y se confunde con las aceleraciones). Se combinan:

El filtro de Kalman extendido (EKF) es la versión sistemática: usa tu modelo para predecir (con el giróscopo y el acelerómetro como entradas) y corrige con GPS, barómetro y magnetómetro, pesando cada fuente con su ruido. PX4 (EKF2) y ArduPilot (EKF3) lo hacen con cuaterniones y error de actitud de 3 componentes, igual que en el capítulo 10. En el curso, la clase 11 (observadores) es la base.

Subir la precisión del modelo

No hace falta un modelo perfecto. Hace falta que sea bueno en lo que afecta al control: alrededor de la frecuencia de corte de los lazos y en los regímenes en los que vas a volar. Ordena las mejoras por lo que cambian en el control, no por lo sofisticadas que parecen.

Los porcentajes son órdenes de magnitud típicos, no garantías.

Medidas baratas que valen mucho

Inercia con péndulo bifilar

Cuelgas el avión de dos hilos paralelos de longitud \(L\) separados \(D\), lo giras un poco alrededor del eje vertical y mides el periodo \(T\) de la oscilación. Repites con el avión colgado en las tres orientaciones.

Banco de empuje

Mide a la vez empuje, par, rpm, corriente y voltaje, barriendo el acelerador. Con eso calibras \(C_T\), \(C_Q\), \(R\), \(I_0\) y \(K_V\) de verdad. Si puedes, repite con viento (túnel o banco en un coche) para tener \(C_T(J)\). Haz también escalones de acelerador para medir \(\tau_m\).

Servos y retrasos

Manda un escalón y un barrido de frecuencia (chirp) a cada servo y mide su posición: sale el ancho de banda y el límite de velocidad. Mide el retraso total del lazo: tiempo entre que el autopiloto lee la IMU y el motor cambia de rpm.

Identificación de sistemas con datos de vuelo

Es la fuente más fiable, pero es un ciclo que se repite, no un último paso:

1. Vuela primero con un controlador prudente, diseñado con el modelo físico.
2. Excita los ejes que te interesan: dobletes, 3-2-1-1 o barridos de frecuencia, en cada mando por separado.
3. Registra a alta frecuencia: IMU, rpm, corriente, posición de los servos, mandos.
4. Estima parámetros físicos (derivadas, \(\tau_m\), retraso) partiendo de tu modelo, no una caja negra. Métodos: error de ecuación (mínimos cuadrados sobre \(\dot p=L_pp+L_{\delta_a}\delta_a+\dots\)), error de salida, o en frecuencia.
5. Valida con otro vuelo distinto del que usaste para identificar.
6. Actualiza el modelo, rediseña, vuelve a comprobar márgenes.

Herramientas: System Identification Toolbox de MATLAB (tfest, ssest, greyest para modelos con parámetros físicos), CIFER (en frecuencia, estándar en helicópteros), o en Python scipy.optimize.least_squares.

Orden recomendado

1. Medir masa, CG, inercias, actuadores y empuje (días de trabajo, gran ganancia).
2. Modelo con strip theory y primeros controladores en simulación.
3. Tablas de VSPAERO y datos de propulsión.
4. CFD o túnel solo donde haga falta (\(\alpha\) alto, interacción hélice-ala).
5. Identificación en vuelo, actualizar el modelo, repetir.

Herramientas: MATLAB, Simulink, Python y simuladores

MATLAB y sus toolboxes

Comprueba lo que tienes con ver. Todo lo de este cuaderno se puede hacer solo con MATLAB base, Control System y Optimization.

Simulink: cómo organizar el modelo

Los mismos bloques del capítulo 1, cada uno en un subsistema: Aerodinámica, Propulsión, Cuerpo rígido (6DOF), Sensores, Estimador, Controlador, Actuadores. Tus funciones .m entran tal cual en bloques MATLAB Function. Ventajas de Simulink: ver las señales, controladores discretos a distinto ritmo, y generar código para el hardware. Inconveniente: esconde las ecuaciones. Recomendación: primero escríbelo en código (para entenderlo), después pásalo a Simulink si te hace falta.

Python

Simuladores y autopilotos

• Gazebo (probablemente lo que oíste como "Gasivo"): simulador de robótica con física y 3D. Se usa con PX4 y ArduPilot en modo SITL (el código real del autopiloto corriendo contra el simulador).
• PX4: autopiloto abierto con soporte para VTOL tail-sitter. Su arquitectura (cascada, mixer, EKF2) es la del capítulo 11.
• ArduPilot: otro autopiloto abierto; los tail-sitters van dentro de su modo QuadPlane.
• JSBSim: el motor de dinámica de vuelo que usan FlightGear y otros. Puedes meter ahí tus tablas.
• FlightGear: visualización 3D gratuita; MATLAB y Simulink le pueden mandar el estado del avión.
• X-Plane: simulador comercial, también se conecta a autopilotos.

Ruta típica hacia el hardware: SIL (controlador como software contra tu modelo) → HIL (el autopiloto real, con su procesador, contra el modelo en el PC) → vuelo atado → vuelo libre.

Herramientas aerodinámicas

• OpenVSP / VSPAERO: geometría y VLM/paneles. Tienes el modelo hecho.
• XFLR5: más sencillo; da derivadas de estabilidad de un avión completo.
• AVL: VLM clásico de Mark Drela, rápido para derivadas.
• XFOIL: perfiles 2D con viscosidad: \(c_{l\alpha}\), \(c_{d0}\), \(c_{l,\max}\).
• CFD: STAR-CCM+ (comercial, lo tienes en la universidad), OpenFOAM y SU2 (abiertos).

Visualización 3D, de menos a más

1. patch + hgtransform en MATLAB con tu geometría: simple y tuyo.
2. Aero.Animation del Aerospace Toolbox.
3. FlightGear conectado a MATLAB o Simulink: escenario realista.
4. Gazebo con PX4 SITL: cuando el autopiloto real entre en juego.

Enlaces y vídeos

• OpenVSP · Base de datos de hélices UIUC
• Documentación de PX4 · ArduPilot · Gazebo
• JSBSim · FlightGear
• python-control · CasADi
• Caughey, Dynamical Equations for Flight Vehicles (referencia de tus clases FD02)
• Vídeos: busca en YouTube "MATLAB Tech Talks Brian Douglas" (series sobre filtros de Kalman, control en la práctica y simulación de drones) y "PX4 tailsitter".

Tu plan de trabajo

Cómo usar la IA sin dejar de aprender

1. Deriva tú primero. Pide a la IA que revise tu derivación, no que la haga.
2. Un archivo cada vez. No aceptes cinco archivos nuevos de golpe. Pide uno y que te explique cada línea con la ecuación que implementa.
3. Predice antes de ejecutar. Antes de correr una simulación, escribe qué esperas. Si no coincide, ahí hay algo que aprender (o un error).
4. Cambia un parámetro y explícalo. Sube la inercia, mueve el CG, quita la estela. Si sabes explicar el cambio, entiendes el modelo.
5. Escribe tú las pruebas de cordura (capítulo 8). Son la mejor forma de saber si el código hace lo que crees.
6. Lleva un diario de decisiones: qué supones, por qué, de dónde sale cada número.

Fases y tareas

Símbolos y referencias

Tus clases, por tema

• Ejes y notación: FD01. Ecuaciones del movimiento: FD02–FD05. Strip theory: FD06 (dos partes).
• Modos del avión: FD13. Cualidades de vuelo: FD14. Ensayos en vuelo y aumento de estabilidad: FD15–FD18.
• Control: modelado (1), estabilidad (3), lugar de las raíces (5), Bode (6), PID (7), espacio de estados (9), diseño de controladores (10), observadores (11), control óptimo (12).

Libros

• M. V. Cook, Flight Dynamics Principles (el libro del curso).
• B. L. Stevens, F. L. Lewis, E. N. Johnson, Aircraft Control and Simulation: el modelo completo con cuaterniones y los coeficientes \(\Gamma\).
• R. W. Beard, T. W. McLain, Small Unmanned Aircraft: Theory and Practice: el más cercano a tu proyecto, de la física al autopiloto.
• B. Etkin, Dynamics of Atmospheric Flight.
• J. G. Leishman, Principles of Helicopter Aerodynamics: cantidad de movimiento y BEMT a fondo.
• M. Drela, Flight Vehicle Aerodynamics.
• V. Klein, E. A. Morelli, Aircraft System Identification; M. B. Tischler, R. K. Remple, Aircraft and Rotorcraft System Identification.