Sistema Integrado de Percepción Visual y Control Adaptativo para Navegación Autónoma de Quadcopters en Entornos Dinámicos Jose Daniel Quintana Fuentes Universidad Tecnológica de Pereira Maestría en Ingeniería Eléctrica Facultad de Ingenierías Pereira, Risaralda, Colombia
2026
Sistema Integrado de Percepción Visual y Control Adaptativo para Navegación Autónoma de Quadcopters en Entornos Dinámicos Jose Daniel Quintana Fuentes Proyecto de grado presentado como requisito parcial para obtener el título de: M.Sc. en Ingeniería Eléctrica Director:
Ph.D. Eduardo Giraldo Suárez Codirector:
MSc. Sergio Velarde Gómez Línea de investigación:
Control Adaptativo, Percepción Visual Grupo de investigación en Control Automático Universidad Tecnológica de Pereira Maestría en Ingeniería Eléctrica Facultad de Ingenierías Pereira, Risaralda, Colombia
2026
Resumen La navegación autónoma de cuadricópteros en entornos tridimensionales dinámicos exige sistemas capaces de percibir, anticipar y evitar obstáculos bajo restricciones de cómputo y actuación. Se plantea un marco de trabajo que integra un modelo dinámico del vehículo, la identificación de un modelo discreto multivariable y un controlador óptimo con acción integral y manejo explícito de saturaciones, complementado con un esquema de adaptación en línea supervisado. Sobre esta base de control se organiza una arquitectura jerárquica de navegación sustentada en mapas tridimensionales de ocupación voxelizada, donde se combinan una planificación global en 3D guiada por muestreo aleatorio, el suavizado geométrico de trayectorias, un planificador local reactivo y un seguidor de trayectorias diseñado para operar en plataformas embebidas. La implementación y evaluación en escenarios simulados con obstáculos estáticos y móviles muestran un comportamiento robusto del sistema de navegación, con trayectorias suaves, libres de colisión y coherentes con las restricciones temporales de la misión. La integración de control óptimo adaptativo y percepción voxelizada ofrece así una base técnica sólida y reproducible para llevar la navegación autónoma de cuadricópteros hacia entornos más complejos y su posterior validación en plataformas físicas.
III
Agradecimientos Quiero expresar mi gratitud a mi familia: mis padres, mi abuela y mi hermana. Gracias por su apoyo incondicional, el cual ha sido la base para poder dar este gran paso. A mi director, gracias por haberme guiado con su paciencia y vasto conocimiento, permitiéndome alcanzar el éxito en este trabajo.
Por último, agradezco al Grupo de Investigación de Control Automático por acogerme con mucho cariño y brindarme un excelente ambiente de trabajo. V
Contenido
1. Introducción
1
1.1.
Planteamiento del problema
. . . . . . . . . . . . . . . . . . . . . . . .
1
1.2.
Justificación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.2.1.
Pertinencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.2.2.
Viabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.2.3.
Impacto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.3.
Objetivos
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.3.1.
Objetivo general
. . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.3.2.
Objetivos específicos . . . . . . . . . . . . . . . . . . . . . . . .
7
1.4.
Estado del arte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2. Modelado Matemático del Cuadricóptero
11
2.1.
Sistemas de Referencia y Convenciones . . . . . . . . . . . . . . . . . .
11
2.2.
Marcos de Referencia y Ángulos de Euler . . . . . . . . . . . . . . . . .
12
2.3.
Formulación Newton–Euler . . . . . . . . . . . . . . . . . . . . . . . . .
14
2.3.1.
Dinámica Traslacional . . . . . . . . . . . . . . . . . . . . . . .
14
2.3.2.
Dinámica Rotacional . . . . . . . . . . . . . . . . . . . . . . . .
15
2.4.
Ecuaciones No Lineales . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
2.5.
Redefinición de las Entradas para Diseño de Control
. . . . . . . . . . .
18
2.6.
Modelo Linealizado . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
2.6.1.
Formalismo general de linealización . . . . . . . . . . . . . . . .
22
2.6.2.
Punto de Operación: Vuelo Estacionario (Hover) . . . . . . . . .
22
2.6.3.
Linealización alrededor de Hover
. . . . . . . . . . . . . . . . .
23
2.7.
Modelo Linealizado en Espacio de Estados
. . . . . . . . . . . . . . . .
24
3. Identificación de Sistemas Multivariables
27
3.1.
Ecuación en diferencias multivariable
. . . . . . . . . . . . . . . . . . .
28
3.2.
Estimación fuera de línea por Mínimos Cuadrados (Identificación Offline)
29
VII
CONTENIDO
3.3.
Estimación en línea . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
3.3.1.
Mínimos Cuadrados Recursivos (RLS)
. . . . . . . . . . . . . .
32
3.4.
Representación de estado extendida
. . . . . . . . . . . . . . . . . . . .
35
3.4.1.
Control en espacio de estados discreto . . . . . . . . . . . . . . .
37
4. Diseño y Análisis del Controlador
39
4.1.
Objetivos de control y restricciones . . . . . . . . . . . . . . . . . . . . .
39
4.2.
Diseño nominal en tiempo discreto . . . . . . . . . . . . . . . . . . . . .
41
4.2.1.
Regulador lineal cuadrático discreto (DLQR) . . . . . . . . . . .
41
4.2.2.
Seguimiento con acción integral . . . . . . . . . . . . . . . . . .
42
4.2.3.
Saturaciones y anti–windup
. . . . . . . . . . . . . . . . . . . .
42
4.3.
Control nominal ante incertidumbre y enlace adaptativo . . . . . . . . . .
43
5. Control Adaptativo Óptimo: STR Indirecto
45
5.1.
Marco General del STR Indirecto . . . . . . . . . . . . . . . . . . . . . .
47
5.2.
Sistema de Identificación . . . . . . . . . . . . . . . . . . . . . . . . . .
49
5.3.
Reajuste discreto del controlador sobre el modelo actualizado . . . . . . .
50
5.4.
Gestión adaptativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
6. Percepción visual para navegación autónoma
53
6.1.
Arquitectura general de percepción y navegación basada en mapas de ocupación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
6.2.
Entorno de simulación y flujo de ejecución . . . . . . . . . . . . . . . . .
55
6.3.
Representación voxelizada del entorno . . . . . . . . . . . . . . . . . . .
56
6.3.1.
Definición y construcción de la grilla de ocupación . . . . . . . .
56
6.3.2.
Mapas voxel precomputados . . . . . . . . . . . . . . . . . . . .
57
6.4.
Planificación global sobre mapas de ocupación tridimensionales
. . . . .
58
6.4.1.
Nodo interactivo de planificación RRT . . . . . . . . . . . . . . .
58
6.4.2.
Formulación geométrica del RRT
. . . . . . . . . . . . . . . . .
59
6.4.3.
Configuración mediante parámetros en YAML
. . . . . . . . . .
60
6.5.
Suavizado de trayectorias mediante B-splines . . . . . . . . . . . . . . .
60
6.5.1.
Nodo de suavizado y algoritmo de Chaikin
. . . . . . . . . . . .
61
6.6.
Seguimiento de trayectorias y máquina de estados de misión
. . . . . . .
62
6.6.1.
Descripción general del seguidor de trayectorias
. . . . . . . . .
62
6.6.2.
Máquina de estados de misión . . . . . . . . . . . . . . . . . . .
63
6.6.3.
Actualización en línea de trayectorias . . . . . . . . . . . . . . .
63
6.7.
Planificación local reactiva y costura de trayectorias . . . . . . . . . . . .
65
VIII
6.7.1.
Ventana local y detección de conflicto . . . . . . . . . . . . . . .
65
6.7.2.
Replanificación local acotada y costura de trayectorias . . . . . .
66
7. Análisis y resultados
69
7.1.
Resultados del modelo ARX MIMO (orden seis, retardo unitario) . . . . .
69
7.1.1.
Métricas de evaluación y validación . . . . . . . . . . . . . . . .
71
7.2.
Resultados de la estimación en línea mediante RLS . . . . . . . . . . . .
72
7.3.
Resultados del controlador nominal LQI . . . . . . . . . . . . . . . . . .
74
7.3.1.
Estructura del lazo y variables empleadas . . . . . . . . . . . . .
75
7.3.2.
Métricas de evaluación y validación . . . . . . . . . . . . . . . .
76
7.4.
Resultados del Control Adaptativo Óptimo — STR Indirecto . . . . . . .
79
7.4.1.
Estructura del esquema STR y variables empleadas . . . . . . . .
79
7.4.2.
Métricas de evaluación y validación . . . . . . . . . . . . . . . .
79
7.4.3.
Referencias aplicadas . . . . . . . . . . . . . . . . . . . . . . . .
80
7.4.4.
Seguimiento de posición y guiñada
. . . . . . . . . . . . . . . .
81
7.4.5.
Arranque seguro con LQI 12×12
. . . . . . . . . . . . . . . . .
81
7.4.6.
Evolución completa de las seis salidas . . . . . . . . . . . . . . .
82
7.4.7.
Errores de seguimiento . . . . . . . . . . . . . . . . . . . . . . .
82
7.4.8.
Lectura conjunta de métricas y supervisión
. . . . . . . . . . . .
82
7.4.9.
Comparación global con controladores de referencia . . . . . . .
84
7.5.
Resultados de navegación y planificación autónoma . . . . . . . . . . . .
86
8. Conclusiones y trabajos futuros
91
Apéndices
92
A. Identificación por Subespacios (MOESP/N4SID)
93
A.1. Fundamento Teórico
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
93
A.2. Implementación y Resultados . . . . . . . . . . . . . . . . . . . . . . . .
94
A.3. Análisis Cuantitativo de Ajuste . . . . . . . . . . . . . . . . . . . . . . .
97
B. Diseño Espectral y Modal por Realimentación de Estados
99
B.1. Asignación de polos en forma canónica
. . . . . . . . . . . . . . . . . . 100
B.2. Asignación de estructuras propias (eigenestructura) . . . . . . . . . . . . 100 B.3. Chequeos estructurales del modelo discreto
. . . . . . . . . . . . . . . . 101
B.3.1. Controlabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 B.3.2. Observabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 IX
CONTENIDO
B.3.3. Resumen numérico para el diseño LQI . . . . . . . . . . . . . . . 103 B.4. Forma canónica de Jordan
. . . . . . . . . . . . . . . . . . . . . . . . . 103
C. Regulador por Realimentación de Estados
109
D. Observadores de Estado (Luenberger / orden reducido)
115
D.1. Diseño del observador Luenberger (orden completo)
. . . . . . . . . . . 115
D.2. Diseño del observador de orden reducido . . . . . . . . . . . . . . . . . . 116 E. Estimación Óptima de Estados: LQE / Filtro de Kalman
121
F. Controladores de Seguimiento de Referencia
125
F.1.
Seguimiento con ganancia de referencia Kr
. . . . . . . . . . . . . . . . 125
F.2.
Seguimiento con acción integral
. . . . . . . . . . . . . . . . . . . . . . 126
G. Sintonía Óptima de Pesos LQR con Algoritmo Genético
129
G.1. Modelo de optimización
. . . . . . . . . . . . . . . . . . . . . . . . . . 130
G.1.1. Diagrama de flujo . . . . . . . . . . . . . . . . . . . . . . . . . . 130 H. Control Lineal Cuadrático Gaussiano (LQG)
135
H.1. Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 Bibliografía
138
X
Figuras
2.1.
Marco inercial {X, Y, Z} con eje Z vertical y vector gravedad −g; el eje X′ ilustra una rotación ψ en el plano horizontal. . . . . . . . . . . . . . .
13
2.2.
Convenciones de ejes y ángulos Euler ZYX: ψ (yaw), θ (pitch) y ϕ (roll). El empuje total actúa a lo largo del eje zb, apuntando hacia arriba en el marco de cuerpo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
2.3.
Plataforma en configuración “X”: numeración de rotores, brazos de longitud l y sentidos de giro. CW indica giro horario y CCW giro antihorario.
17
3.1.
Señales de excitación PRBS en configuración MIMO por canal de entrada.
34
3.2.
Diagrama de bloques del esquema de identificación RLS–ARX MIMO en Simulink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
4.1.
Esquema general del controlador LQI. a) representación compacta del lazo integral con el regulador LQR; b) realización detallada en tiempo discreto en espacio de estados. . . . . . . . . . . . . . . . . . . . . . . . . .
43
5.1.
Esquema conceptual del flujo Modelado →Identificación →Control Óptimo →Control Adaptativo (STR indirecto). Las flechas sólidas representan el flujo nominal entre capítulos; las flechas discontinuas ilustran el bucle de adaptación con datos en línea y rediseño periódico de ganancias.
47
5.2.
Arquitectura general de un regulador autoajustable (STR) indirecto. El módulo de identificación de sistema estima los parámetros del proceso y el bloque de diseño actualiza en línea los parámetros del controlador que actúa sobre la planta. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
XI
FIGURAS
6.1.
Arquitectura general del subsistema de percepción y navegación basada en mapas de ocupación, implementada en ROS Melodic sobre una NVI- DIA Jetson Nano. El mapa voxelizado se comparte entre los módulos de planificación global, suavizado de trayectorias, seguimiento y planificación local. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
6.2.
Vista combinada del entorno de simulación en Gazebo y su representación voxelizada en RViz. A la derecha se observa la disposición física de los obstáculos, mientras que a la izquierda se muestra el mapa 3D discretizado utilizado por los módulos de planificación y navegación. . . . . . .
56
6.3.
Discretización de una nube de puntos tridimensional en una malla de vóxeles con resolución fija. Cada punto se proyecta sobre una celda de la grilla, que se marca como ocupada si pertenece al dominio del mapa. . . .
57
6.4.
Ejemplo de ruta global generada por el planificador RRT sobre el mapa de ocupación voxelizado del laboratorio. La trayectoria poligonal en rojo conecta el punto inicial con el objetivo evitando las regiones ocupadas representadas por los bloques azules. . . . . . . . . . . . . . . . . . . . .
59
6.5.
Ruta poligonal generada por el planificador RRT y trayectoria suavizada obtenida a partir del algoritmo de Chaikin en el mismo entorno.
. . . . .
62
6.6.
Máquina de estados del seguidor de trayectorias. Se representan los estados idle, taking_off, following y landing, junto con las transiciones asociadas a la recepción de nuevas trayectorias y a las condiciones de vuelo.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
64
6.7.
Esquema de planificación local en una ventana tridimensional alrededor del cuadricóptero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
6.8.
Secuencia de replanificación local y costura de trayectorias: (i) ruta global suavizada, (ii) detección del tramo conflictivo en la ventana local, (iii) generación de un camino local alternativo mediante RRT restringido y (iv) trayectoria resultante tras la costura. . . . . . . . . . . . . . . . . . .
67
7.1.
Superposición medida–predicha (predicción de un paso) para el modelo ARX MIMO de orden seis con retardo unitario, estimado por mínimos cuadrados fuera de línea. . . . . . . . . . . . . . . . . . . . . . . . . . .
70
7.2.
Estimación en línea (RLS) para ARX MIMO de orden seis y retardo unitario — predicción a un paso. . . . . . . . . . . . . . . . . . . . . . . . .
72
7.3.
Errores absolutos de estimación por salida en el esquema RLS–ARX MIMO. 73
XII
7.4.
Errores relativos de estimación por salida ( %) en el esquema RLS–ARX
MIMO.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
7.5.
Esquema resumido del lazo LQI nominal utilizado en la evaluación de resultados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
7.6.
Seguimiento de referencias en x, y, z, ψ con LQI. . . . . . . . . . . . . .
76
7.7.
Entradas de control u = [T, τx, τy, τz]⊤durante el seguimiento. . . . . . .
77
7.8.
Errores de seguimiento ex, ey, ez, eψ. . . . . . . . . . . . . . . . . . . . .
78
7.9.
Esquema resumido del STR indirecto utilizado en la evaluación de resultados.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
80
7.10. Referencias seleccionadas para x, y, z y ψ (multi–step). . . . . . . . . . .
82
7.11. Seguimiento de x, y, z (arriba) y ψ (abajo). Se sobreponen referencias
(líneas punteadas).
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
7.12. Arranque seguro con LQI 12×12: comparación salida–referencia por canal. 84
7.13. Evolución de las 6 salidas del cuadricóptero (posición y actitud). . . . . .
85
7.14. Errores ex, ey, ez y eψ a lo largo de la prueba. . . . . . . . . . . . . . . .
86
7.15. Secuencia cronológica de navegación autónoma en un entorno interior di-
námico. (i) Detección de conflicto sobre la trayectoria activa y solicitud de replaneo local. (ii) Recepción de una trayectoria alternativa para evasión del obstáculo. (iii) Seguimiento de la trayectoria actualizada tras el replaneo. (iv) Continuación de la misión hacia el objetivo de navegación. .
87
A.1. Comparación entre datos medidos y predichos utilizando el método N4SID. 95 A.2. Comparación entre datos medidos y predichos utilizando el método MOESP. 96 B.1. Comparación de polos en lazo cerrado para el cuadricóptero: ambas técnicas alcanzan el mismo espectro deseado. . . . . . . . . . . . . . . . . . 101 C.1. Comparación de respuesta al escalón (todas las salidas x, y, z, ϕ, θ, ψ). Arriba: lazo abierto. Abajo: lazo cerrado con u = −Kx.
. . . . . . . . . 111
C.2. Respuesta al escalón con realimentación (u = [1, 1, 1, 1]T). Salidas de posición: x, y, z.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
C.3. Respuesta al escalón con realimentación (u = [1, 1, 1, 1]T). Salidas angulares: ϕ, θ, ψ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 C.4. Salidas x, y, z, ϕ, θ, ψ superpuestas ante u = [1, 1, 1, 1]T. La línea punteada marca una referencia de comparación fija, utilizada sólo como guía visual. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
XIII
FIGURAS
D.1. Observador de Luenberger (orden completo). Comparación x vs. ˆx para los estados 1–6. El transitorio breve confirma la estabilidad de A −LC. . 116 D.2. Observador de Luenberger (orden completo). Comparación x vs. ˆx para los estados 7–12. La inyección de innovaciones corrige las variables no medidas directamente.
. . . . . . . . . . . . . . . . . . . . . . . . . . . 117
D.3. Observador de orden reducido. Salidas medidas y (línea continua) frente a Cˆx (línea discontinua). La coincidencia respalda la ubicación de polos del estimador. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 E.1. LQE: comparación x vs. ˆx para los estados 1–6, (x, y, z, ˙x, ˙y, ˙z). Se observa un seguimiento cercano con transitorios suaves aun en presencia de ruido. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 E.2. LQE: comparación x vs. ˆx para los estados 7–12, (ϕ, θ, ψ, ˙ϕ, ˙θ, ˙ψ). La corrección por innovación alinea las pendientes y reduce el sesgo residual. 123 F.1.
Seguimiento multi–step con ganancia de referencia Kr: todas las salidas {x, y, z, ϕ, θ, ψ} frente a referencias a trozos constantes. La referencia se inyecta simultáneamente en los seis canales; Kr corrige la ganancia de DC y K gobierna el transitorio. . . . . . . . . . . . . . . . . . . . . . . . 126 F.2.
Seguimiento multi–step con acción integral: todas las salidas y sus referencias a trozos. El integrador elimina error estacionario ante cambios de consigna y sesgos persistentes. . . . . . . . . . . . . . . . . . . . . . . . 127 G.1. Flujo del AG para sintonía de Q, R en LQR. . . . . . . . . . . . . . . . . 131 G.2. LQR sintonizado por AG (multi–step): (x, y, z) vs referencia. . . . . . . . 133 G.3. LQR sintonizado por AG (multi–step): (ϕ, θ, ψ) vs referencia.
. . . . . . 134
H.1. LQG: escalón unitario en seis salidas. Transitorios acotados en (x, y, z) y variables angulares próximas a cero por diseño del costo. . . . . . . . . . 136 H.2. LQG: seguimiento multi–step. Los ángulos se conservan cercanos a cero y sólo intervienen transitoriamente para sostener la traslación.
. . . . . . 137
XIV
Capítulo 1 Introducción
1.1.
Planteamiento del problema En los últimos años, el mercado de los quadcopters ha experimentado un crecimiento significativo a nivel mundial y en Colombia. Según un informe reciente de Marketsand- Markets, se espera que el mercado global de UAVs (Vehículos Aéreos No Tripulados) crezca de 30.2 mil millones de dólares en 2024 a 48.5 mil millones de dólares para 2029, con una tasa de crecimiento anual compuesta (CAGR) del 9.9 % [1]. Además, se proyecta que el mercado global de drones alcanzará un valor de 54.6 mil millones de dólares en 2030, con el mercado comercial creciendo a una tasa de crecimiento anual compuesta (CAGR) del 7.7 % [2].
Por otra parte, la NASA ha proyectado una expansión continua en la industria de UAVs, con aplicaciones que van desde la entrega de paquetes y la agricultura de precisión hasta misiones de exploración en Marte [3]. Este crecimiento es impulsado por la creciente demanda de soluciones tecnológicas innovadoras y sostenibles en diversos sectores industriales y comerciales a nivel global. Dichas innovaciones están alineadas con los Objetivos de Desarrollo Sostenible (ODS) de las Naciones Unidas, especialmente el ODS 9, que promueve la construcción de infraestructuras resilientes, la industrialización inclusiva y sostenible, y el fomento de la innovación [4].
La Federación Internacional de Robótica (IFR) ha destacado el creciente uso de drones en América Latina, y específicamente en Colombia, debido a su potencial para optimizar la eficiencia operativa en diversos sectores económicos [5]. En línea con esta tendencia, la implementación de estas aeronaves no tripuladas en el país se ha expandido rápidamente en áreas como la agricultura, la seguridad y la logística. Según un informe de la Cámara Colombiana de Comercio Electrónico (CCCE), el mercado de drones en el país
CAPÍTULO 1. INTRODUCCIÓN
ha mostrado un crecimiento anual del 22 % [6], lo cual se encuentra alineado con el Plan Nacional de Desarrollo 2022–2026 del gobierno colombiano, que prioriza la adopción de tecnologías emergentes para fortalecer la productividad y la sostenibilidad económica [7]. Además, Statista proyecta que el mercado de drones en Colombia alcanzará un valor de
8.5 millones de euros en 2024, con una tasa de crecimiento anual del 6.0 % entre 2024 y
2028; asimismo, se espera que el volumen del mercado llegue a 23.3 mil unidades para 2028, con un incremento del 7.0 % en 2025 [8].
Este panorama refleja un escenario en el que la demanda de drones capaces de operar con mayor autonomía es cada vez más relevante. A medida que estas aeronaves se involucran en aplicaciones industriales más complejas y exigentes, como la inspección de infraestructura o la logística de última milla, resulta imperativo mejorar sus capacidades de navegación y evasión de obstáculos en entornos dinámicos. Esta mejora no solo consolidaría su valor tecnológico y económico, sino que también contribuiría a optimizar la eficiencia operativa en los sectores que impulsan su creciente adopción [9]. Considerando el creciente interés internacional y las proyecciones favorables en mercados como el colombiano, resulta pertinente enmarcar el rol de estas aeronaves en el ecosistema tecnológico actual. En la era de las Industrias 4.0, la automatización y la robótica están transformando el panorama a una velocidad acelerada. En este escenario, los vehículos autónomos emergen como una innovación revolucionaria, con el potencial de impactar sectores tan diversos como la agricultura, la medicina, la defensa y la exploración espacial [10, 11, 12, 13]. Esta tecnología, en constante evolución, se clasifica según el medio de desplazamiento, abarcando vehículos autónomos terrestres (AV), acuáticos (AUV) y aéreos (UAV) [14]. Dentro de estos últimos, la principal ventaja es su capacidad para alcanzar rápidamente ubicaciones deseadas, superando las limitaciones temporales y geográficas propias de los vehículos terrestres y acuáticos [10]. Los UAV, comúnmente denominados drones, pueden dividirse adicionalmente según su aerodinámica y características físicas [15, 16, 17]. Entre estas clasificaciones, los quadcópteros, un tipo de multirrotor, destacan por su estabilidad, versatilidad y facilidad de control. Gracias a estas cualidades, se han consolidado como una opción preferente en aplicaciones que van desde la topografía y la inspección de infraestructura hasta la búsqueda y rescate en espacios reducidos, así como la entrega de mercancías [15]. Al conjugar la creciente demanda con las capacidades intrínsecas de estas plataformas aéreas, se refuerza la relevancia de optimizar su autonomía y su habilidad para navegar de manera segura en entornos complejos.
Para cumplir con las funciones de las aplicaciones descritas, es necesario garantizar
la seguridad e integridad del dron a través de la navegación autónoma, de forma que el vehículo pueda desplazarse por un entorno sin intervención humana [18, 10]. Además, se debe realizar una planificación de rutas que le permita desenvolverse en entornos desconocidos, donde no se cuenta con información previa sobre los posibles obstáculos ni sus dinámicas de movimiento [12, 19]. Esta planificación de rutas es conocida como "Local Path Planning", y su implementación eficiente es un problema latente debido a su importancia y complejidad [20, 12, 21]. Para lograrlo, es fundamental que el vehículo tenga la capacidad de percibir su entorno en tiempo real mientras navega [22, 23]. Los sensores son las herramientas a través de las cuales el dron percibe su entorno. Por este motivo, la elección adecuada de sensores es una cuestión clave para los UAV: su limitada capacidad de carga impone restricciones significativas debido al peso adicional que estos implican [18, 24, 15, 10], lo que puede resultar en una reducción de la duración de vuelo y una menor capacidad de respuesta ante cambios en el entorno [25]. La clasificación de los sensores se determina según el método a través del cual obtienen la información del entorno [10, 26]. Aquellos que emiten una señal propia para la detección o extracción de información se denominan sensores activos, e incluyen dispositivos tales como el LiDAR, Sonar y Radar [10, 11]. Estos sensores cuentan con la ventaja de no verse afectados por la iluminación; sin embargo, presentan limitaciones al momento de detectar objetos transparentes o de pequeño tamaño, además de que su precisión a larga distancia se ve afectada [24]. Por el contrario, aquellos dispositivos que capturan la información del entorno utilizando la percepción visual, como las cámaras convencionales, de profundidad y de eventos, son sensores pasivos [27, 24]. Si bien estos sensores presentan un menor consumo de energía respecto a los sensores activos y ofrecen una detección detallada, son sensibles a condiciones de luz y clima y pueden presentar limitaciones con la tasa de captura, generando latencia según la cámara utilizada, lo que limita su velocidad de captura de movimiento [28, 10, 27].
Los métodos de percepción visual emplean cámaras y técnicas de visión por computadora para identificar y analizar elementos presentes en el entorno, con el propósito de mejorar la navegación autónoma en vehículos aéreos no tripulados (UAVs) [29, 24, 30, 31]. Estas técnicas posibilitan el reconocimiento y seguimiento de objetos, la estimación de sus posiciones y trayectorias, así como la integración con sensores adicionales para optimizar la obtención de información. Diversas estrategias han sido exploradas para operar en entornos desconocidos, incluyendo aquellas que implementan enfoques tradicionales de procesamiento de imagen, algoritmos de extracción de características y métodos basados en aprendizaje profundo [32, 33, 34, 35, 36, 37]. No obstante, la alta variabilidad
CAPÍTULO 1. INTRODUCCIÓN
en las condiciones ambientales y la necesidad de respuestas en tiempo real plantean retos significativos en cuanto a robustez, eficiencia y precisión, lo que ha motivado la búsqueda de herramientas más especializadas.
En este contexto, la detección y clasificación de objetos mediante Redes Neuronales Convolucionales (CNN) ofrece niveles de precisión sobresalientes en escenarios complejos [30, 31]; sin embargo, su implementación suele requerir datos especializados, mayor tiempo de entrenamiento y una cuidadosa sincronización de fuentes. La segmentación de imágenes a través de arquitecturas especializadas, como U-Net o Mask R-CNN, permite identificar con detalle regiones de interés y distinguir con nitidez los contornos de los objetos [28, 38]; no obstante, estas arquitecturas incrementan la demanda computacional y pueden dificultar la ejecución sobre plataformas aéreas con recursos limitados [32, 35]. Por su parte, el seguimiento de objetos (Object Tracking), al incorporar modelos temporales como las Redes Neuronales Recurrentes (RNN, por sus siglas en inglés Recurrent Neural Networks) o arquitecturas de memoria a corto y largo plazo (LSTM, por sus siglas en inglés Long Short-Term Memory), posibilita estimar poses y predecir trayectorias con mayor fidelidad [29, 24]; sin embargo, esta integración exige un balance delicado entre el procesamiento espacio-temporal y la necesidad de mantener la eficiencia en tiempo real. Estas aproximaciones resultan prometedoras para entornos dinámicos y heterogéneos; no obstante, la optimización de sus componentes y el diseño de arquitecturas livianas siguen siendo retos latentes que motivan la exploración de nuevas soluciones [37]. A pesar de que las mejoras en percepción visual y control autónomo han incrementado la capacidad de reacción ante estímulos inmediatos del entorno [39, 40, 41], la ausencia de una perspectiva anticipatoria limita el desempeño en entornos dinámicos y altamente inciertos. Muchos de los sistemas actuales se basan en un enfoque predominantemente reactivo, centrándose en las condiciones presentes sin contemplar posibles cambios futuros. Esta dependencia de la información instantánea introduce demoras en la toma de decisiones y reduce la precisión en las maniobras, particularmente cuando los obstáculos alteran con rapidez su posición o presentan trayectorias complejas. No obstante, en [42, 43, 44] se proponen métodos que integran capacidades anticipatorias y adaptativas, lo que sugiere un camino hacia la consolidación de enfoques donde la detección, la predicción del movimiento de obstáculos y el control adaptativo converjan para mejorar la eficiencia y la resiliencia de la navegación autónoma en tiempo real. A su vez, los planteamientos tradicionales de control, tales como los métodos clásicos y los enfoques de control óptimo, entre ellos el regulador cuadrático lineal (LQR, por sus siglas en inglés Linear Quadratic Regulator), han mostrado su eficacia en escena-
rios donde las condiciones se mantienen relativamente estables y bien definidas [34], sin embargo, la falta de mecanismos de identificación y estimación en tiempo real implica que estos enfoques resulten menos efectivos cuando el entorno presenta una alta variabilidad y la dinámica del sistema se torna incierta [40]. En contraste, las estrategias de control adaptativo y predictivo ofrecen un margen de maniobra más amplio, ajustando parámetros según el contexto y fortaleciendo la resiliencia del dron ante alteraciones súbitas [45, 46]. Esta transición hacia esquemas más flexibles y anticipatorios no sólo responde a la necesidad de enfrentar condiciones cambiantes, sino que también se alinea con el panorama previamente delineado, donde la percepción y la predicción del entorno convergen con el control, abriendo paso a sistemas de navegación autónoma más integrales y versátiles. Considerando estas circunstancias, surge la siguiente pregunta de investigación: ¿Cómo desarrollar un sistema integrado de percepción visual y control adaptativo para mejorar la adaptabilidad y seguridad en la navegación autónoma de un quadcopter para la evasión de obstáculos en entornos dinámicos?
1.2.
Justificación
1.2.1.
Pertinencia En la era de las Industrias 4.0, la automatización y la robótica continúan revolucionando múltiples sectores, desde la agricultura hasta la logística, la seguridad y la exploración espacial [10, 11, 12, 13]. En este contexto, los vehículos aéreos no tripulados (UAVs), y en particular los quadcopters, han demostrado un notable potencial gracias a su versatilidad, estabilidad y capacidad para desenvolverse con precisión en entornos complejos [10]. Sin embargo, el reto de asegurar una navegación verdaderamente autónoma y eficiente persiste, especialmente en escenarios donde la dinámica de los obstáculos se vuelve impredecible, modificando sus trayectorias de manera rápida y constante [12, 19]. Esta necesidad de anticipación y ajuste continuo ante condiciones cambiantes orienta el desarrollo de un sistema capaz de integrar percepción visual y control adaptativo, combinando la visión por computador con enfoques de control predictivo. Con ello, se busca superar las limitaciones propias de sistemas reactivos convencionales, incorporando estrategias que no solo reconozcan y sigan los movimientos de los obstáculos, sino que también predigan su comportamiento futuro [37, 43, 42]. Esta aproximación, soportada en técnicas como el controlador adaptativo LQR complementado con un filtro de Kalman, aspira a proporcionar a los drones la habilidad de ajustar sus maniobras en tiempo real, incrementando así su eficiencia y seguridad en entornos dinámicos [26, 45].
CAPÍTULO 1. INTRODUCCIÓN
Al consolidar una navegación más robusta, este enfoque no solo amplía su aplicabilidad en misiones críticas —tales como el rescate en zonas de difícil acceso, el monitoreo de áreas urbanas complejas o la optimización de procesos agrícolas avanzados—, sino que también impulsa un progreso tangible en el desarrollo científico-técnico de la navegación autónoma de UAVs. Esta evolución no se limita al ámbito tecnológico, pues al integrarse en sectores estratégicos y responder a demandas emergentes, contribuye a la construcción de entornos más resilientes y sostenibles. En este sentido, se alinea con objetivos globales como el ODS 9 de la ONU, orientado hacia la innovación y la industrialización sostenible [4]. De este modo, las soluciones resultantes no solo atienden desafíos inmediatos, sino que refuerzan la posición de las tecnologías autónomas en un mercado dinámico y en constante transformación, sentando las bases para un impacto perdurable.
1.2.2.
Viabilidad El proyecto se fundamenta en la vasta experiencia y el sólido conocimiento del grupo de investigación en el campo de los vehículos aéreos no tripulados (UAVs). Este tema está a la vanguardia de la tecnología, lo que ha llevado a numerosos grupos de investigación en todo el mundo a dedicar grandes esfuerzos para desarrollar soluciones innovadoras en la navegación autónoma de UAVs [47]. Un aspecto crucial del desarrollo de este proyecto es el respaldo de avances significativos logrados por el grupo de investigación en control automático de la Universidad Tecnológica de Pereira, bajo la dirección del profesor Eduardo Giraldo Suárez. Este grupo ha trabajado en diversas investigaciones y desarrollos, como el diseño e implementación de controladores multivariable en un UAV cuadricóptero para fumigación aérea [48], el desarrollo de control robusto acoplado de sistemas multivariables aplicado a UAVs [49], y el control multivariable no lineal aplicado a UAVs [50]. Además, han investigado la implementación de técnicas de control inteligente en un helicóptero no tripulado de dos grados de libertad [51]. Estas investigaciones demuestran un historial consolidado en el diseño y control de UAVs, asegurando que el grupo posee el conocimiento y las herramientas necesarias para llevar a cabo con éxito el presente proyecto.
1.2.3.
Impacto Con el desarrollo de este proyecto de investigación, se busca avanzar en el estudio y mejora de la seguridad y eficiencia en la navegación autónoma de UAVs, lo que reducirá el riesgo de colisiones y aumentará la capacidad de respuesta ante obstáculos imprevistos.
Este tipo de estudios tiene un gran impacto en diversas aplicaciones relacionadas en el área del desarrollo tecnológico [47]. Además, ampliará las aplicaciones de los drones en áreas como la monitorización ambiental, inspección de infraestructuras y misiones de búsqueda y rescate [15]. Este avance contribuirá al progreso tecnológico en el ámbito de los vehículos autónomos, proporcionando nuevas metodologías y enfoques que pueden ser adoptados por otros investigadores y desarrolladores.
1.3.
Objetivos
1.3.1.
Objetivo general Desarrollar un sistema integrado de percepción visual y control adaptativo para mejorar la adaptabilidad y seguridad en la navegación autónoma de un quadcopter para la evasión de obstáculos en entornos dinámicos.
1.3.2.
Objetivos específicos Desarrollar un sistema de percepción visual para la detección temprana de obstáculos dinámicos.
Desarrollar una estrategia de control adaptativo que integre los datos del sistema de percepción visual para ajustar las maniobras del quadcopter.
Evaluar el desempeño del sistema integrado mediante simulaciones y pruebas experimentales en entornos dinámicos controlados.
1.4.
Estado del arte En las últimas décadas, la investigación en vehículos aéreos no tripulados (UAVs) ha generado una amplia variedad de enfoques y soluciones orientadas a incrementar la autonomía, adaptabilidad y seguridad de estas plataformas. Este interés se debe, en gran medida, a la necesidad de operar en entornos complejos y cambiantes, donde los UAVs deben reaccionar eficientemente ante obstáculos dinámicos y condiciones ambientales adversas [32, 33, 34, 35, 36, 41, 52, 19]. Entre las plataformas más estudiadas se encuentran los quadcopters, cuya versatilidad y estabilidad los han posicionado como casos de estudio idóneos para la exploración de técnicas de percepción visual, predicción del entorno y control adaptativo.
CAPÍTULO 1. INTRODUCCIÓN
En el campo de la percepción, la visión por computadora ha sido esencial para la detección de objetos, clasificación de imágenes y segmentación de escenas, empleando sensores pasivos como cámaras RGB, de profundidad o de eventos [27, 28]. Las Redes Neuronales Convolucionales (CNN) han demostrado una eficacia notable en la identificación y clasificación de obstáculos, apoyando la navegación en entornos urbanos densamente poblados [30, 29], mientras que estrategias basadas en flujo óptico, tales como el algoritmo Lucas-Kanade, han ofrecido robustez ante movimientos de alta frecuencia [27, 26]. No obstante, estos enfoques enfrentan limitaciones en cuanto a requerimientos computacionales, dependencia de condiciones ambientales favorables y dificultades para responder en tiempo real sobre plataformas aéreas con recursos limitados [32, 35]. Además de los enfoques orientados a detección, segmentación y seguimiento, en la literatura se ha consolidado una línea de investigación centrada en la integración directa entre percepción visual y control, conocida como control servovisual (visual servoing) [53, 54, 55]. En estos métodos, la información visual extraída de la escena se emplea dentro del lazo de control para guiar el movimiento del sistema, ya sea a partir de características definidas en el plano imagen o de variables reconstruidas en el espacio tridimensional [53, 54]. Esta línea ha sido ampliamente estudiada en robótica y vehículos autónomos, y ha dado lugar a aplicaciones en tareas de seguimiento, guiado, alineación y aterrizaje autónomo [56, 57, 58]. No obstante, su implementación también enfrenta retos asociados a la calibración, al campo de visión, a la sensibilidad ante perturbaciones visuales y a las exigencias de procesamiento en tiempo real [55, 58].
Aparte de la detección instantánea, la capacidad de anticipar dinámicas futuras se ha convertido en un factor determinante para la operación en entornos dinámicos. En este sentido, la incorporación de modelos temporales como las Redes Neuronales Recurrentes (RNN) o arquitecturas de memoria a corto y largo plazo (LSTM, por sus siglas en inglés Long Short-Term Memory) ha permitido el seguimiento de trayectorias, el modelado temporal y la predicción del movimiento de obstáculos [29, 24]. Esta capacidad anticipatoria, al proveer información sobre las posibles trayectorias futuras, resulta especialmente valiosa cuando se integra con estrategias de control adaptativo y predictivo, ya que proporciona un insumo crítico para ajustar las maniobras antes de que las alteraciones ocurran [39, 40, 42].
En cuanto al control, las estrategias clásicas, como los controladores PID o el Regulador Cuadrático Lineal (LQR, por sus siglas en inglés Linear Quadratic Regulator), han sido ampliamente empleadas gracias a su relativa sencillez de implementación y su efectividad en condiciones operativas bien definidas. Por ejemplo, el LQR, especialmente
cuando se combina con técnicas de estimación como el Filtro de Kalman, proporciona un equilibrio entre eficiencia y precisión al regular el comportamiento del dron en entornos con variabilidad moderada [34, 43]. Sin embargo, ante las limitaciones en entornos con alta incertidumbre, donde las dinámicas no se ajustan a modelos lineales sencillos, el interés se ha volcado hacia estrategias más flexibles. En estos casos, la dependencia de modelos estáticos y la ausencia de mecanismos internos de ajuste limitan el desempeño de los enfoques clásicos [40].
En contraste, las estrategias de control adaptativo y predictivo han emergido como alternativas capaces de superar estas deficiencias. Estas metodologías incorporan información en tiempo real sobre el estado del dron, las características del entorno y las posibles trayectorias de los obstáculos para ajustar los parámetros de control de forma continua y anticipada [38, 44, 45, 46]. De esta manera, no solo corrigen desviaciones en la trayectoria, sino que también anticipan alteraciones futuras, lo que les confiere una mayor capacidad para enfrentar condiciones cambiantes y escenarios dinámicos. Además, la integración con datos de percepción visual y predicción del entorno fortalece la sinergia entre detección y maniobra, facilitando la convergencia hacia sistemas más resilientes y robustos ante perturbaciones externas [42, 43, 41].
La implementación y validación de estas soluciones ha venido acompañada de metodologías de evaluación en simuladores de alta fidelidad, entornos controlados y esquemas de Hardware-in-the-Loop (HIL), en los cuales componentes físicos interactúan con entornos virtuales simulados, permitiendo probar el desempeño del sistema bajo condiciones cercanas a la realidad sin incurrir en riesgos significativos [59, 60, 61, 62, 63, 64]. Estas prácticas facilitan el análisis de la robustez ante condiciones adversas, la eficiencia de los algoritmos de percepción y predicción, así como la precisión del control autónomo. Dichas limitaciones ponen de manifiesto la necesidad de sistemas integrados capaces de unir percepción, predicción y control adaptativo en un solo marco metodológico, abordando así el desafío central planteado en el presente trabajo. De esta forma, el estado del arte revela el crecimiento de un cuerpo de conocimiento que combina percepción visual avanzada, predicción dinámica y control adaptativo, evidenciando las oportunidades y retos pendientes para consolidar sistemas integrales capaces de navegar en tiempo real con mayor autonomía, adaptabilidad y seguridad.
La Tabla 1.1 sintetiza las principales líneas de investigación identificadas en la literatura, junto con las metodologías empleadas y sus limitaciones.
CAPÍTULO 1. INTRODUCCIÓN
Tabla 1.1. Resumen del estado del arte. Línea Referentes principales y enfoques Limitaciones Percepción visual
• [27, 28]: cámaras RGB, de profundidad
y de eventos.
• [30, 29]: detección y clasificación
mediante CNN.
• [26]: flujo óptico y seguimiento visual.
• Alta demanda computacional.
• Sensibilidad a iluminación y
condiciones ambientales.
• Dificultad para operación en tiempo
real sobre plataformas embebidas. Control servovisual
• [53, 54, 55]: fundamentos de visual
servoing basados en imagen y en pose.
• [56, 57, 58]: aplicaciones recientes en
guiado, seguimiento y aterrizaje autónomo.
• Dependencia de calibración y
geometría de visión.
• Restricciones de campo de visión.
• Sensibilidad a perturbaciones visuales
y exigencias de procesamiento en tiempo real.
Predicción del entorno
• [29, 24]: modelado temporal mediante
RNN y LSTM.
• [39, 40, 42]: predicción de trayectorias
y anticipación de maniobras.
• Mayor complejidad computacional.
• Compromiso entre precisión predictiva
y tiempo de respuesta.
• Integración no siempre directa con el
control.
Control clásico y óptimo
• [34, 43]: controladores PID, LQR e
integración con filtro de Kalman.
• Regulación efectiva en condiciones
nominales y moderadamente variables.
• Dependencia de modelos fijos.
• Menor robustez ante incertidumbre y
cambios dinámicos.
• Limitada capacidad de autoajuste.
Control adaptativo y predictivo
• [38, 44, 45, 46]: control adaptativo,
predictivo y ajuste continuo de parámetros.
• [42]: integración entre percepción,
predicción y maniobra.
• Mayor complejidad de implementación.
• Requiere integración robusta entre
percepción, estimación y control.
• Validación experimental más exigente.
Validación experimental
• [59, 60, 61]: esquemas HIL y
validación en tiempo real.
• [62, 63, 64]: simuladores de alta
fidelidad y plataformas de prueba para UAVs.
• La simulación no reproduce toda la
complejidad del entorno real.
• Persisten brechas entre validación
virtual y despliegue físico.
Capítulo 2 Modelado Matemático del Cuadricóptero Se desarrolla el modelado matemático no lineal de un cuadricóptero (quadcopter) de 6 grados de libertad (6DOF) con acoplamientos giroscópicos. El enfoque se basa en la formulación Newton–Euler para derivar las ecuaciones de movimiento tanto traslacionales como rotacionales. Se muestra el procedimiento algebraico que conduce a las ecuaciones finales, y se organiza el modelo en forma de espacio de estados no lineal. Este modelo habilita el desarrollo de estrategias de control basadas en la dinámica no lineal y, adicionalmente, permite obtener un modelo lineal alrededor de la condición de hover, entendida como el vuelo estacionario en el que el vehículo mantiene aproximadamente constantes su posición y orientación, para facilitar el diseño y análisis de controladores lineales.
2.1.
Sistemas de Referencia y Convenciones El cuadricóptero posee cuatro rotores dispuestos en forma de “X”. Cada rotor genera un empuje fi y, debido a la inercia de los rotores, se producen acoplamientos giroscópicos. Se definen las siguientes variables:
x, y, z: Posición del centro de masa en el marco inercial {I}. ˙x, ˙y, ˙z: Velocidades lineales.
ϕ, θ, ψ: Ángulos de Euler (roll, pitch, yaw).
˙ϕ, ˙θ, ˙ψ: Velocidades angulares.
fi: Empuje generado por el rotor i (con fi = k ω2 i ).
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
ωi: Velocidad angular del rotor i.
m: Masa del cuadricóptero.
g: Aceleración gravitatoria.
Ixx, Iyy, Izz: Momentos de inercia del vehículo. IR: Momento de inercia de cada rotor (para acoplamientos giroscópicos). l: Brazo (distancia desde el centro de masa a cada rotor).
2.2.
Marcos de Referencia y Ángulos de Euler Marcos.
En este modelado se adopta un sistema de ejes ortonormal, es decir, X×Y = Z y xb × yb = zb:
{X, Y, Z} (marco inercial): se asume Z apuntando hacia arriba, de modo que la gravedad actúa en la dirección [0, 0, −g]⊤, es decir, en −Z, como se muestra en la Figura 2.1.
{xb, yb, zb} (marco de cuerpo): se adhiere al cuadricóptero y se define con ejes:
• xb: hacia el frente del vehículo,
• yb: hacia la izquierda,
• zb: hacia arriba.
Con esta elección, el producto cruz xb × yb apunta en la dirección +zb, de acuerdo con la regla de la mano derecha. Para el marco inercial, se ha definido Z hacia arriba, de modo que la fuerza gravitatoria actúa en la dirección −Z. Orientación (Euler ZYX).
Para describir la orientación, se utilizan tres ángulos de Euler:
ϕ (roll), θ (pitch), ψ (yaw), donde:
ϕ describe la rotación alrededor del eje xb. θ describe la rotación alrededor del eje yb. ψ describe la rotación alrededor del eje zb.
X Y Z −g Marco Inercial {X, Y, Z} ψ X′ Figura 2.1. Marco inercial {X, Y, Z} con eje Z vertical y vector gravedad −g; el eje X′ ilustra una rotación ψ en el plano horizontal. Adoptando la convención Euler ZYX, la matriz de rotación que transforma vectores del marco de cuerpo {B} al marco inercial {I} se escribe como: R(ϕ, θ, ψ) = Rz(ψ) Ry(θ) Rx(ϕ), con vI = RvB y, por ortonormalidad, R−1 = R⊤. En esta expresión, las rotaciones elementales están dadas por:
Rx(ϕ) =
1
0
0
0
cos ϕ −sen ϕ
0
sen ϕ cos ϕ
,
Ry(θ) = cos θ
0
sen θ
0
1
0
−sen θ
0
cos θ
,
Rz(ψ) = cos ψ −sen ψ
0
sen ψ cos ψ
0
0
0
1
.
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
De este modo, la matriz de rotación completa resulta:
R(ϕ, θ, ψ) = cos ψ cos θ sen ϕ sen θ cos ψ −cos ϕ sen ψ sen ϕ sen ψ + cos ϕ sen θ cos ψ sen ψ cos θ sen ϕ sen ψ sen θ + cos ϕ cos ψ −sen ϕ cos ψ + cos ϕ sen ψ sen θ −sen θ sen ϕ cos θ cos ϕ cos θ
.
En esta parametrización, la orientación del cuerpo se obtiene mediante una rotación de yaw ψ alrededor del eje Z inercial, seguida de un pitch θ alrededor del eje Y ′ y, finalmente, un roll ϕ alrededor del eje X′′. La dependencia trigonométrica en ϕ, θ, ψ introduce la no linealidad en la transformación y determina cómo se proyectan los vectores (por ejemplo, la fuerza de empuje) del marco de cuerpo {B} al marco inercial {I}. Como se ilustra en la Figura 2.2, cada rotor genera un empuje individual fi dirigido a lo largo del eje zb; la contribución conjunta de estos empujes da lugar al empuje total del vehículo. Figura 2.2. Convenciones de ejes y ángulos Euler ZYX: ψ (yaw), θ (pitch) y ϕ (roll). El empuje total actúa a lo largo del eje zb, apuntando hacia arriba en el marco de cuerpo.
2.3.
Formulación Newton–Euler
2.3.1.
Dinámica Traslacional Fuerza de Empuje en el Marco de Cuerpo Sea fi el empuje de cada rotor. La suma de empujes se define como T = f1 + f2 + f3 + f4.
Dado que en {B} (marco de cuerpo) se ha definido zb apuntando hacia arriba, la fuerza de empuje total actúa a lo largo de +zb. En coordenadas del marco de cuerpo, el vector de
empuje se expresa como Fb =
0
0
T , físicamente el empuje neto apunta en la misma dirección que el eje zb. Transformación al Marco Inercial Para llevar Fb al marco inercial {I}, se multiplica por la matriz de rotación R(ϕ, θ, ψ) que transforma vectores de {B} a {I}. Además, se resta la fuerza de la gravedad (que actúa en −Z con Z hacia arriba en el marco inercial): FI = R(ϕ, θ, ψ) Fb −
0
0
m g .
La fuerza resultante FI se aplica en F = m a para obtener las aceleraciones lineales ¨x, ¨y, ¨z en el marco inercial.
Componentes de la Aceleración Definiendo la posición del centro de masa en el marco inercial como p = [x, y, z]T, se cumple:
m ¨p = FI.
La multiplicación R(ϕ, θ, ψ) Fb involucra las funciones trigonométricas de ϕ, θ, ψ. Cada componente (X, Y, Z) de la fuerza transformada aporta un término a ¨x, ¨y, ¨z.
2.3.2.
Dinámica Rotacional Ecuación de Euler La dinámica rotacional del cuadricóptero se describe mediante la ecuación de Euler:
I ˙ωB + ωB × I ωB = τ, donde ωB = [p, q, r]T es el vector de velocidades angulares expresado en el marco de cuerpo {B}, τ es el vector de torques aplicados y I es la matriz de inercia del cuadricóptero respecto a su centro de masa. Se adopta la forma diagonal I = diag(Ixx, Iyy, Izz), donde Ixx, Iyy e Izz son los momentos de inercia respecto a los ejes xb, yb y zb, respectivamente.
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
Torques en Configuración en “X” Para un cuadricóptero en configuración “X”, tal como se ilustra en la Figura 2.3, los rotores giran en sentidos alternados: CW (clockwise, sentido horario) y CCW (counterclockwise, sentido antihorario), vistos desde arriba. Esta disposición permite compensar, en condiciones nominales, el par de arrastre aerodinámico total y establecer el torque de guiñada a partir del desequilibrio entre rotores contrarrotantes. Bajo esta configuración, los torques de roll (τϕ), pitch (τθ) y yaw (τψ) se definen combinando diferencialmente los empujes:
τϕ = l√
2
f2 + f3 −f1 −f4
,
τθ = l√
2
f1 + f2 −f3 −f4
,
τψ = b ω2
1 −ω2
2 + ω2
3 −ω2
4
,
y se puede añadir un término IR ˙Ωpara los efectos giroscópicos. Relación entre (p, q, r) y ( ˙ϕ, ˙θ, ˙ψ) La matriz que relaciona las velocidades angulares del marco de cuerpo con las derivadas de Euler es: p q r =
1
0
−sen θ
0
cos ϕ sen ϕ cos θ
0
−sen ϕ cos ϕ cos θ ˙ϕ ˙θ ˙ψ .
Despejando ˙p, ˙q, ˙r en la ecuación de Euler y sustituyendo τϕ, τθ, τψ, se obtienen las ecuaciones para ¨ϕ, ¨θ, ¨ψ.
2.4.
Ecuaciones No Lineales Dinámica traslacional Sea T = f1 + f2 + f3 + f4 la fuerza total generada por los rotores. En el sistema de cuerpo {B}, cuya orientación se define con xb hacia adelante, yb hacia la izquierda y zb hacia arriba, el vector de empuje se expresa como Fb = [ 0, 0, T]T. Al transformarlo al sistema inercial {I}, con eje Z positivo hacia arriba, y considerando que la gravedad actúa en −Z, se obtienen las siguientes ecuaciones de movimiento traslacional:
¨x = T m sen ϕ sen ψ + sen θ cos ϕ cos ψ
,
(2.1)
¨y = T m −sen ϕ cos ψ + sen ψ sen θ cos ϕ
,
(2.2)
¨z = −g + T m cos ϕ cos θ.
(2.3)
Front
1
2
3
4
CCW
CCW
CW CW l Figura 2.3. Plataforma en configuración “X”: numeración de rotores, brazos de longitud l y sentidos de giro. CW indica giro horario y CCW giro antihorario. El término −g aparece porque en el sistema inercial {I} el eje Z se orienta hacia arriba, mientras que la gravedad actúa en sentido opuesto. Los términos proporcionales a T/m corresponden a la proyección del eje zb, donde actúa el empuje en el marco de cuerpo, sobre {I}, lo que coincide con la tercera columna de R(ϕ, θ, ψ). Dinámica rotacional Para la dinámica rotacional se adopta una configuración en “X”, en la cual el marco de cuerpo {B} se define con xb orientado al frente, yb hacia la izquierda y zb hacia arriba, y los brazos se encuentran a ±45◦respecto a xb. En este caso, los torques se modelan como: τϕ = l√
2
f2 + f3 −f1 −f4
,
τθ = l√
2
f1 + f2 −f3 −f4
,
τψ = b ω2
1 −ω2
2 + ω2
3 −ω2
4
+ término giroscópico
,
donde l es la semidistancia al centro, b el coeficiente de par aerodinámico (arrastre rotacional), ωi las velocidades angulares de los rotores y τgyro el posible efecto giroscópico debido a la inercia de los mismos. Como fi ∝ω2 i , los torques dependen directamente de las velocidades de giro de los motores.
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
Bajo el supuesto de un tensor de inercia diagonal I = diag(Ixx, Iyy, Izz), la dinámica rotacional del vehículo se describe mediante las ecuaciones de Euler: Ixx ¨ϕ + (Izz −Iyy) ˙θ ˙ψ = τϕ,
(2.4)
Iyy ¨θ + (Ixx −Izz) ˙ϕ ˙ψ = τθ,
(2.5)
Izz ¨ψ + (Iyy −Ixx) ˙ϕ ˙θ = τψ.
(2.6)
Cinemática de Euler El vínculo entre las velocidades angulares del cuerpo ωb = [ p, q, r ]T y las tasas de cambio de los ángulos de Euler η = [ ϕ, θ, ψ ]T se establece mediante: ˙η = W(ϕ, θ) ωb, W(ϕ, θ) =
1
sen ϕ tan θ cos ϕ tan θ
0
cos ϕ −sen ϕ
0
sen ϕ/ cos θ cos ϕ/ cos θ .
La relación inversa, útil para simulación, es:
ωb = T(ϕ, θ) ˙η, T(ϕ, θ) =
1
0
−sen θ
0
cos ϕ sen ϕ cos θ
0
−sen ϕ cos ϕ cos θ .
Debe señalarse que existe una singularidad cinemática en cos θ = 0, correspondiente a θ = ±90◦.
2.5.
Redefinición de las Entradas para Diseño de Control En el diseño de control de cuadricópteros es habitual reagrupar los empujes individuales {f1, f2, f3, f4} en entradas virtuales {u1, u2, u3, u4}. Esta redefinición simplifica la síntesis de control, pues cada entrada se asocia de forma directa a un grado de libertad específico del vehículo. Adoptando la convención de marco de cuerpo: xb: hacia el frente, yb: hacia la izquierda, zb: hacia arriba,
se define:
u1 = f1 + f2 + f3 + f4, u2 = l √
2
f2 + f3 −f1 −f4
,
u3 = l √
2
f1 + f2 −f3 −f4
,
u4 = b ω2
1 −ω2
2 + ω2
3 −ω2
4
o, de forma equivalente, u4 = τψ.
Interpretación:
u1: fuerza total de empuje, que interviene en la dinámica traslacional (u1 m en ¨x, ¨y, ¨z).
u2 y u3: torques netos de roll y pitch, obtenidos a partir de combinaciones diferenciales de los empujes.
u4: torque neto de yaw, derivado en la diferencia de momentos aerodinámicos entre rotores.
Esta agrupación es especialmente útil porque permite asignar un único comando a cada grado de libertad: u1 controla la altitud, u2 y u3 la inclinación (roll y pitch), y u4 la orientación en yaw.
Dinámica Traslacional Sea T = u1 la fuerza total de empuje. Al proyectar esta fuerza sobre el marco inercial mediante la matriz de rotación R(ϕ, θ, ψ), se obtiene:
¨x = u1 m sen ϕ sen ψ + sen θ cos ϕ cos ψ
,
(2.7)
¨y = u1 m −sen ϕ cos ψ + sen ψ sen θ cos ϕ
,
(2.8)
¨z = −g + u1 m cos ϕ cos θ.
(2.9)
Interpretación:
¨x: aceleración en la dirección de avance, modulada por la inclinación del dron. ¨y: aceleración lateral, dependiente de la orientación en roll y yaw. ¨z: aceleración vertical, que combina el efecto de la gravedad con la proyección del empuje sobre el eje Z inercial.
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
Dinámica Rotacional Aplicando la formulación de Newton–Euler y las redefiniciones de torques, la dinámica rotacional queda descrita por:
Ixx ¨ϕ + (Izz −Iyy) ˙θ ˙ψ = u2,
(2.10)
Iyy ¨θ + (Ixx −Izz) ˙ϕ ˙ψ = u3,
(2.11)
Izz ¨ψ + (Iyy −Ixx) ˙ϕ ˙θ = u4 + IR ˙Ω,
(2.12)
donde:
u2 y u3: torques netos en roll y pitch.
u4: torque neto de yaw, expresado como τψ o a partir de las velocidades angulares de los rotores.
IR ˙Ω: efecto giroscópico asociado a la inercia de los rotores y la velocidad neta de giro Ω.
Relación entre Velocidades Angulares y Ángulos de Euler La cinemática de Euler establece el vínculo entre las velocidades angulares del cuerpo ωb = (p, q, r)T y las derivadas de los ángulos de Euler ( ˙ϕ, ˙θ, ˙ψ)T. Este mapeo se expresa mediante:
p q r =
1
0
−sen θ
0
cos ϕ sen ϕ cos θ
0
−sen ϕ cos ϕ cos θ ˙ϕ ˙θ ˙ψ , donde se observa que las velocidades angulares en el marco de cuerpo son combinaciones lineales de las tasas de cambio de los ángulos de Euler, con dependencia no lineal en ϕ y θ. Cabe señalar que esta parametrización presenta una singularidad en θ = ±90◦. Formulación No Lineal para el Diseño de Control Con la redefinición de entradas virtuales (u1, u2, u3, u4), las ecuaciones de movimiento del cuadricóptero pueden reescribirse de forma compacta de la siguiente manera:
Dinámica traslacional:
¨x = u1 m sen ϕ sen ψ + sen θ cos ϕ cos ψ
,
(2.13)
¨y = u1 m −sen ϕ cos ψ + sen ψ sen θ cos ϕ
,
(2.14)
¨z = −g + u1 m cos ϕ cos θ.
(2.15)
Dinámica rotacional:
Ixx ¨ϕ + (Izz −Iyy) ˙θ ˙ψ = u2,
(2.16)
Iyy ¨θ + (Ixx −Izz) ˙ϕ ˙ψ = u3,
(2.17)
Izz ¨ψ + (Iyy −Ixx) ˙ϕ ˙θ = u4 + IR ˙Ω.
(2.18)
La introducción de las entradas virtuales permite asignar de forma directa cada ui a un grado de libertad particular:
u1: controla la altitud, al representar la fuerza total de empuje proyectada sobre los ejes inerciales.
u2 y u3: ajustan la inclinación en roll y pitch, respectivamente. u4: regula el yaw, asociado a la rotación alrededor del eje vertical. Esta formulación es especialmente ventajosa para la síntesis de controladores, ya que desacopla conceptualmente cada entrada en su grado de libertad correspondiente y simplifica el diseño de leyes de control para el cuadricóptero.
2.6.
Modelo Linealizado Se considera un modelo lineal aproximado del cuadricóptero, obtenido a partir de la linealización del modelo no lineal, con el fin de facilitar el diseño de controladores clásicos y modernos. La linealización se realiza alrededor del punto de operación correspondiente al vuelo estacionario (hover), empleando una expansión de Taylor de primer orden y aproximaciones de ángulo pequeño.
La convención utilizada es consistente con las secciones previas: {I} marco inercial con Z ↑, {B} marco de cuerpo con xb al frente, yb hacia la izquierda, zb ↑; parametrización de ángulos de Euler ZYX (roll ϕ, pitch θ, yaw ψ); y entradas virtuales {u1, u2, u3, u4}.
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
2.6.1.
Formalismo general de linealización Sea el sistema de múltiples entradas y múltiples salidas (MIMO, por sus siglas en inglés Multiple-Input Multiple-Output) no lineal en espacio de estados: ˙x = f(x, u), y = h(x, u),
(2.19)
con estado x ∈Rn, entrada u ∈Rm, salida y ∈Rr. Sea (x0, u0) un punto de equilibrio tal que f(x0, u0) = 0. Definiendo perturbaciones δx = x−x0, δu = u−u0, la expansión de Taylor de primer orden alrededor del equilibrio es: ˙δx = ∂f ∂x
(x0,u0) | {z } A δx + ∂f ∂u
(x0,u0) | {z } B δu + O(∥δx, δu∥2),
(2.20)
δy = ∂h ∂x
(x0,u0) | {z } C δx + ∂h ∂u
(x0,u0) | {z } D δu + O(∥δx, δu∥2).
(2.21)
Despreciando términos de orden superior, se obtiene el modelo linealizado en espacio de estados:
˙δx = A δx + B δu, δy = C δx + D δu.
(2.22)
2.6.2.
Punto de Operación: Vuelo Estacionario (Hover) En vuelo estacionario, el vehículo mantiene posición y orientación constantes con velocidades nulas. En términos de las variables de estado: ˙x = ˙y = ˙z = 0, ¨x = ¨y = ¨z = 0, ϕ = θ = ψ = 0, ˙ϕ = ˙θ = ˙ψ = 0,
(2.23)
El equilibrio vertical exige que el empuje total compense el peso: u1,0 = T0 = m g, u2,0 = u3,0 = u4,0 = 0.
(2.24)
Si cada rotor produce empuje cuadrático fi = kf ω2 i , entonces en equilibrio ω1,0 = ω2,0 = ω3,0 = ω4,0 = ω0,
4 kf ω2
0 = m g ⇒ω0 =
smg 4kf
.
(2.25)
Las condiciones (2.23)–(2.25) son consistentes con el modelo no lineal: con ϕ = θ = 0 se obtiene ¨z = −g + u1/m, de modo que u1 = mg implica ¨z = 0.
Tabla 2.1. Punto de operación en vuelo estacionario (hover). Variable Valor Comentario ϕ, θ, ψ
0
actitud nula p, q, r
0
tasas nulas ˙x, ˙y, ˙z
0
velocidad nula u1,0 mg empuje total en equilibrio u2,0, u3,0, u4,0
0
torques nulos ω1,0 = · · · = ω4,0 q mg/(4kf) si fi = kfω2 i Tabla 2.2. Parámetros físicos del vehículo (notación y unidades). Símbolo Descripción Unidad m masa kg Ixx, Iyy, Izz momentos principales de inercia kg m2 l semidistancia al centro m kf coef. de empuje (fi = kfω2 i ) N s2 kτ coef. de par (τψ,i = kτω2 i ) N m s2 g aceleración de la gravedad m/s2
2.6.3.
Linealización alrededor de Hover Definiendo perturbaciones δx y δu alrededor de (x0, u0), y aplicando el desarrollo de Taylor de primer orden, se adoptan las siguientes hipótesis de linealización: Aproximaciones de ángulo pequeño:
sen ϕ ≈ϕ, cos ϕ ≈1, sen θ ≈θ, cos θ ≈1, sen ψ ≈ψ, cos ψ ≈1.
Empuje en equilibrio: en el punto de operación se cumple T0 = mg (esto es, u1,0 = mg), mientras que u2,0 = u3,0 = u4,0 = 0. Truncamiento a primer orden: se descartan productos de perturbaciones (p. ej., ϕ θ, θ δu1, ˙ϕ ˙ψ).
En este contexto, y con el fin de resaltar los alcances y limitaciones inherentes al modelo aproximado, resulta conveniente puntualizar el siguiente aspecto:
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
Aproximaciones de ángulo pequeño y validez local Para |ϕ|, |θ|, |ψ| ≪1 (rad), se adoptan: sen α ≃α, cos α ≃1. La precisión del modelo decrece fuera de este rango y cerca de la singularidad cos θ = 0.
Bajo estas premisas, y en concordancia con la cinemática Euler ZYX previamente introducida, la relación exacta ωb = T(ϕ, θ) ˙η (con η = [ϕ θ ψ]T) se evalúa en ϕ = θ = 0; en consecuencia, T(0, 0) = I3 y se obtiene la cinemática linealizada δ ˙ϕ = δp, δ ˙θ = δq, δ ˙ψ = δr.
(2.26)
Por otra parte, al considerar la dinámica traslacional no lineal y proyectar el empuje a {I} mediante la tercera columna de R(ϕ, θ, ψ), la evaluación en hover junto con las aproximaciones de ángulo pequeño conduce, tras retener únicamente términos de primer orden, a δ¨x = g δθ, δ¨y = −g δϕ, δ¨z = 1 m δu1.
(2.27)
Así, las aceleraciones horizontales quedan acopladas linealmente con las actitudes δϕ, δθ: bajo esta convención, una perturbación positiva de pitch δθ > 0 genera una aceleración positiva δ¨x hacia adelante, mientras que una perturbación positiva de roll δϕ > 0 produce una aceleración lateral negativa δ¨y < 0. La aceleración vertical depende directamente de la variación de empuje total.
De modo análogo, partiendo de las ecuaciones de Euler con I = diag(Ixx, Iyy, Izz) y despreciando productos de tasas, se llega a la dinámica rotacional linealizada δ ¨ϕ = 1 Ixx δu2, δ¨θ = 1 Iyy δu3, δ ¨ψ = 1 Izz δu4 + IR Izz δ ˙Ω.
(2.28)
El término adicional IR Izz δ ˙Ωrecoge, en primera aproximación, el efecto giroscópico debido a la inercia de los rotores.
En síntesis, las expresiones (2.26)–(2.28) constituyen el núcleo del modelo lineal local alrededor de hover; a partir de ellas, y mediante la construcción estándar en espacio de estados, se derivan directamente las matrices A, B, C y D empleadas en la representación en espacio de estados que se introduce a continuación.
2.7.
Modelo Linealizado en Espacio de Estados Con el propósito de disponer de una representación idónea para síntesis y análisis, se adopta la descripción en espacio de estados en torno al equilibrio de hover, trabajando en todo momento con variables de desviación respecto del punto de operación. En
consecuencia, se definen x = δx δy δz δ ˙x δ ˙y δ ˙z δϕ δθ δψ δ ˙ϕ δ ˙θ δ ˙ψ T
,
u = h δu1 δu2 δu3 δu4 iT .
Bajo una expansión de Taylor de primer orden alrededor del punto (x0, u0), la dinámica y la salida adoptan la forma ˙x = A x + B u, y = C x + D u, donde A ∈R12×12, B ∈R12×4, C ∈R6×12 y D ∈R6×4. Asimismo, las matrices del modelo se obtienen rigurosamente como los jacobianos de la dinámica y de la función de salida, evaluados en el equilibrio (x0, u0):
A = ∂f(x, u) ∂x
(x0,u0)
,
B = ∂f(x, u) ∂u
(x0,u0)
,
C = ∂h(x, u) ∂x
(x0,u0)
,
D = ∂h(x, u) ∂u
(x0 A efectos de claridad expositiva, se presenta en primer término la forma compacta por bloques, que condensa ceros e identidades sin pérdida de contenido. Con la convención adoptada y tras evaluar los jacobianos en hover, se obtiene:
A =
03
I3
03
03
03
03
0
g
0
−g
0
0
0
0
0
03
03
03
03
I3
03
03
03
03
,
B = 05×4
1
m
0
0
0
0
1
Ixx
0
0
0
0
1
Iyy
0
0
0
0
1
Izz
.
En lo que respecta a la definición de las salidas del sistema, y con el propósito de mantener coherencia con las variables habitualmente medidas por la instrumentación de a bordo, resulta técnicamente adecuado seleccionar y = h δx δy δz δϕ δθ δψ iT , de modo que C actúa como selector de dichas componentes y la salida no depende directamente de la entrada; en consecuencia, C = I3
03
03
03
03
03
I3
03
, D = 06×4.
Esta elección no es la única posible: si se requiriesen velocidades lineales o tasas angulares en la salida, bastaría con redefinir C para incluirlas, manteniendo la coherencia con el punto de operación.
CAPÍTULO 2. MODELADO MATEMÁTICO DEL CUADRICÓPTERO
Finalmente, y únicamente con propósito de referencia exhaustiva, se presenta a continuación la expresión expandida de A y B (sin abreviaciones). Dado que C y D son selectores y nulos, respectivamente, su forma compacta previa resulta suficiente para la discusión siguiente.
A =
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
−g
0
0
0
0
0
0
0
0
0
0
−g
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
,
B =
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
m
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
Ixx
0
0
0
0
1
Iyy
0
0
0
0
1
Izz
.
El modelo linealizado obtenido constituye una aproximación local del comportamiento del cuadricóptero en torno al vuelo estacionario, válida para perturbaciones de pequeña magnitud respecto al punto de equilibrio. En consecuencia, las propiedades de controlabilidad y observabilidad dependen de esta condición de operación y pueden variar de forma significativa en escenarios de vuelo más agresivos o alejados del régimen estacionario.
Capítulo 3 Identificación de Sistemas Multivariables La identificación de sistemas multivariables (MIMO) en este trabajo persigue un objetivo aplicado: construir modelos discretos verificables que sirvan de base directa para la síntesis de control. En particular, se requiere un modelo nominal para el diseño en tiempo discreto (Capítulo 4) y un mecanismo de actualización paramétrica que habilite el reajuste adaptativo durante la operación (Capítulo 5). Con este propósito, el capítulo presenta la formulación autorregresiva con entradas exógenas (ARX, por sus siglas en inglés AutoRegressive with eXogenous inputs) para sistemas MIMO, su estimación fuera de línea mediante Mínimos Cuadrados y su actualización en línea mediante Mínimos Cuadrados Recursivos (RLS, por sus siglas en inglés Recursive Least Squares), dando lugar a una representación en espacio de estados utilizable en el diseño del controlador. El alcance se sitúa en un régimen cercano al vuelo estacionario y muestreo uniforme, con entradas definidas por empuje y momentos y salidas provistas por el sistema de estimación disponible. Se modela la dinámica lineal local y los acoplamientos relevantes; efectos como saturaciones, ruido y pequeñas no linealidades se tratan como perturbaciones. Los detalles operativos de adquisición y preprocesamiento se integran de forma sucinta en las secciones de estimación y se documentan cuantitativamente en el capítulo de resultados.
Sobre esta base, la relación entrada–salida se organiza en forma de regresión lineal ARX, que permite: (i) estimar parámetros en bloque con Mínimos Cuadrados a partir de datos previamente adquiridos; (ii) actualizar dichos parámetros en operación mediante RLS con posibles factores de olvido; y (iii) convertir la regresión a una representación de estado extendida que entrega matrices (Fa, G, Ca) utilizables en el diseño LQR/LQI y en
CAPÍTULO 3. IDENTIFICACIÓN DE SISTEMAS MULTIVARIABLES
la construcción de observadores. Las verificaciones estructurales mínimas (controlabilidad/observabilidad) se aplican como condición de uso del modelo; la evidencia cuantitativa y las métricas de desempeño se presentan más adelante en el capítulo de resultados.
3.1.
Ecuación en diferencias multivariable En esta sección se fija el marco formal para la identificación MIMO en tiempo discreto. Consideramos un sistema con p salidas y q entradas, descrito mediante el operador de retardo r−1 (también denotado q−1 en la literatura). La relación entrada–salida se expresa como A(r−1) y[k] = B(r−1) u[k],
(3.1)
donde y[k] ∈Rp×1 y u[k] ∈Rq×1 son los vectores de salida y entrada respectivamente, y los polinomios matriciales A(·) y B(·) se definen por A(r−1) = A0 + A1r−1 + A2r−2 + · · · + Anr−n,
(3.2)
B(r−1) = B1r−1 + B2r−2 + · · · + Bmr−m,
(3.3)
con n ≥m, Ai ∈Rp×p y Bi ∈Rp×q. En notación explícita, y[k] = y1[k] y2[k]
...
yp[k]
,
u[k] = u1[k] u2[k]
...
uq[k]
.
(3.4)
En la práctica de identificación ARX, se adopta la hipótesis de polinomio A(q) mónico, esto es, A0 = Ip (matriz identidad), lo cual permite reescribir (3.1) como h Ip + A1r−1 + · · · + Anr−ni y[k] = h B1r−1 + · · · + Bmr−mi u[k],
(3.5)
o, en forma recursiva, y[k] = −A1y[k −1] −· · · −Any[k −n] + B1u[k −1] + · · · + Bmu[k −m]. (3.6) Las matrices de coeficientes Ai y Bi recogen los acoplamientos cruzados entre canales:
Ai = ai
11
ai
12
· · ·
ai 1p ai
21
ai
22
· · ·
ai 2p
...
...
...
...
ai p1 ai p2
· · ·
ai pp
,
Bi = bi
11
bi
12
· · ·
bi 1q bi
21
bi
22
· · ·
bi 2q
...
...
...
...
bi p1 bi p2
· · ·
bi pq
.
(3.7)
Así, el i-ésimo componente de la salida puede escribirse, a partir de (3.6)–(3.7), como yi[k] = − n X ℓ=1 p X j=1 aℓ ij yj[k −ℓ] + m X ℓ=1 q X j=1 bℓ ij uj[k −ℓ],
(3.8)
lo que evidencia la dependencia de cada salida respecto a pasados de todas las salidas y entradas (estructura MIMO).
Para fines de estimación, (3.6) puede organizarse en forma de regresión lineal: y[k] = Θ⊤ϕ[k −1], k ≥0,
(3.9)
donde Θ ∈R(pn+qm)×p apila los parámetros de A(·) y B(·), y ϕ[k −1] ∈R(pn+qm)×1 es el vector regresor que reúne salidas y entradas pasadas. En forma compacta, ϕ[k −1] = h y⊤[k −1]
· · ·
y⊤[k −n] u⊤[k −1]
· · ·
u⊤[k −m] i⊤,
(3.10)
y la matriz de parámetros transpuesta adopta la estructura Θ⊤= h −A1
· · ·
−An B1
· · ·
Bm i
,
(3.11)
entendida como concatenación por bloques columna (coherente con el orden de ϕ). Esta forma deja explícito que la identificación ARX se reduce a estimar Θ en un problema lineal en los parámetros.
La representación ARX MIMO y su reorganización en forma de regresión (3.9)–(3.10) constituyen la base operativa para la estimación de parámetros a partir de datos. Los órdenes (n, m) y el orden efectivo por canal se determinarán mediante criterios de información y validación cruzada; con ello se implementarán los esquemas de estimación por Mínimos Cuadrados (fuera de línea) y RLS (en línea), apoyados en un diseño de excitaciones que garantice identificabilidad en el régimen de operación del cuadricóptero. La conversión ARX →espacio de estados proveerá realizaciones discretas utilizables en el diseño nominal del Capítulo 4. En la siguiente sección se desarrollan los fundamentos y las decisiones prácticas que articulan este proceso.
3.2.
Estimación fuera de línea por Mínimos Cuadrados (Identificación Offline) El procedimiento de estimación fuera de línea tiene como propósito obtener los parámetros del modelo multivariable a partir de un conjunto de datos experimentales, empleando el método de Mínimos Cuadrados Ordinarios (OLS). En este enfoque, toda la
CAPÍTULO 3. IDENTIFICACIÓN DE SISTEMAS MULTIVARIABLES
información de entrada y salida se recopila previamente y se utiliza en bloque para calcular los coeficientes del modelo que mejor reproducen la dinámica observada del sistema. Considerando la ecuación en diferencias multivariable descrita en la Subsección 3.1, el modelo puede expresarse como:
y[i] = −A1y[i−1]−· · ·−Any[i−n]+B1u[i−1]+· · ·+Bmu[i−m] = ˆθTϕ[i−1], (3.12) donde ˆθ es el vector de parámetros estimados y ϕ[i −1] es el vector regresor que contiene entradas y salidas pasadas. En forma matricial:
Y = Φˆθ + E,
(3.13)
donde Y agrupa las observaciones de salida, Φ contiene los regresores formados por los datos históricos de entrada y salida, y E representa el vector de errores de predicción. El objetivo es determinar ˆθ que minimice la función de costo cuadrática: V (θ) = 1
2ETE = 1
2(Y −Φˆθ)T(Y −Φˆθ),
(3.14)
cuyo gradiente nulo conduce a la ecuación normal:
ΦTΦˆθ = ΦTY.
(3.15)
Si la matriz ΦTΦ es no singular, la estimación óptima por Mínimos Cuadrados se obtiene como:
ˆθ = (ΦTΦ)−1ΦTY.
(3.16)
Esta expresión corresponde a un esquema de estimación fuera de línea, en el que el cálculo de los parámetros se efectúa a partir del conjunto completo de datos disponible. En este trabajo, la estimación se realizó con las series de entrada y salida previamente adquiridas, obteniéndose el vector de parámetros por mínimos cuadrados de acuerdo con
(3.16).
La metodología se aplicó sobre registros experimentales tomados en condiciones controladas de operación, empleando maniobras excitatorias acotadas que preservan márgenes de seguridad. El modelo resultante se adopta como referencia nominal para el diseño de controladores (Capítulo 4) y como punto de partida para la identificación adaptativa en línea mediante RLS (Capítulo 5).
Esquema procedimental A partir de la formulación ARX multivariable establecida en 3.1 y su organización como regresión lineal en (3.9)–(3.10), así como del estimador de Mínimos Cuadrados
Algorithm 1 Estimación fuera de línea (OLS) para un modelo ARX MIMO de orden seis Require: u ∈RN×4, y ∈RN×6; órdenes (na, nb, nk) = (6, 6, 1) 1: k0 ←m´ax(na, nb + nk) + 1; Nv ←N −k0 + 1 2: Construir Φ ∈RNv×(6na+4nb) concatenando, para k = k0:N:
[−y(k −1), . . . , −y(k −na)] y [ u(k −nk), . . . , u(k −nk −nb + 1)] 3: Definir Y ←y(k0:N) 4: Estimar ˆθ resolviendo m´ınθ∥Y −Φθ∥2 2:
(ecuaciones normales) ˆθ = (Φ⊤Φ)−1Φ⊤Y (estable numéricamente) ˆθ = Φ†Y vía QR/SVD (forma regularizada) (Φ⊤Φ + εI)−1Φ⊤Y 5: Predicción a un paso: ˆY ←Φˆθ; residuo: E ←Y −ˆY (3.12)–(3.16), conviene condensar el flujo operativo que enlaza datos, construcción del regresor, estimación y verificación inicial. De este modo, la ecuación en diferencias se traduce en un procedimiento ordenado que, primero, prepara y alinea el conjunto de muestras; después, forma el regresor Φ apilando retardos de salidas y entradas; a continuación, estima los parámetros en bloque; y, finalmente, valida mediante predicción de un paso, preparando el terreno para los contrastes cuantitativos posteriores. En este marco, el modelo adoptado es un ARX multivariable de orden seis con retardo unitario: cada componente de la salida actual depende linealmente de seis retardos de todas las salidas y de seis retardos de todas las entradas, mientras que el primer retardo en la rama exógena fija el retardo unitario ((3.6)–(3.10)). La elección de órdenes determina el índice inicial k0 y el subconjunto válido de muestras, tras lo cual la estimación por mínimos cuadrados, eventualmente regularizada para mejorar el acondicionamiento, proporciona ˆθ y habilita la predicción ˆY = Φˆθ. Por tanto, la comparación directa entre Y y ˆY complementada con métricas agregadas y el examen cualitativo del residuo ofrece una lectura transparente de la fidelidad del ajuste, mientras que la reagrupación de ˆθ en bloques { ˆAℓ, ˆBℓ} facilita el análisis estructural y su posterior uso tanto en la conversión a espacio de estados como en la actualización recursiva mediante RLS.
3.3.
Estimación en línea Mientras que la estimación fuera de línea requiere disponer de todo el conjunto de datos antes de calcular los parámetros, en muchos sistemas reales especialmente aquellos que operan en tiempo real los modelos deben actualizarse continuamente conforme se
CAPÍTULO 3. IDENTIFICACIÓN DE SISTEMAS MULTIVARIABLES
reciben nuevas mediciones. Este enfoque da lugar a los esquemas de estimación en línea (online), donde los parámetros del modelo se ajustan de manera iterativa a partir de la información disponible en cada instante de muestreo.
En este contexto, el vector de parámetros estimados ˆθ[k] en el instante k se actualiza a partir de su valor anterior ˆθ[k −1] mediante una ley de corrección proporcional al error de estimación:
ˆθ[k] = ˆθ[k −1] + M[k −1] e[k],
(3.17)
donde M[k −1] es la matriz de ganancia del algoritmo y e[k] representa el error de estimación definido como e[k] = y[k] −ˆy[k],
(3.18)
siendo ˆy[k] = ˆθ[k −1]Tϕ[k −1] la salida estimada a partir del vector regresor ϕ[k −1]. El principio esencial de estos métodos radica en actualizar los parámetros del modelo de manera adaptativa, de modo que la estimación refleje cambios graduales en la dinámica del sistema. Este mecanismo resulta especialmente útil en el caso del cuadricóptero, donde las variaciones en las condiciones de vuelo, el desgaste de los actuadores o perturbaciones externas pueden modificar los parámetros efectivos del modelo. La estimación en línea permite así mantener una representación actualizada de la planta y soportar estrategias de control adaptativo, como las desarrolladas en el Capítulo 5.
A continuación, se presenta el algoritmo de Mínimos Cuadrados Recursivos (RLS), uno de los esquemas de estimación en línea más robustos y ampliamente utilizados.
3.3.1.
Mínimos Cuadrados Recursivos (RLS) El algoritmo de Mínimos Cuadrados Recursivos surge de extender el método clásico de Mínimos Cuadrados hacia una formulación recursiva, donde los parámetros se actualizan en cada paso a medida que nuevos datos se incorporan. Esta aproximación evita el procesamiento en bloque de todo el conjunto de datos y reduce significativamente el coste computacional, manteniendo la precisión del estimador.
Partiendo de la expresión general de la estimación fuera de línea (3.16), se define una forma recursiva en la cual las matrices involucradas se actualizan incrementalmente. Sea P[k] ∈R(pn+qm)×(pn+qm) la matriz de covarianza asociada a la incertidumbre de los parámetros, la actualización recursiva de los parámetros se expresa como: ˆθ[k] = ˆθ[k −1] + P[k −1]ϕ[k −1]
1 + ϕ[k −1]TP[k −1]ϕ[k −1]
y[k] −ϕ[k −1]T ˆθ[k −1]
, (3.19)
y la actualización de la matriz de covarianza:
P[k] = P[k −1] −P[k −1]ϕ[k −1]ϕ[k −1]TP[k −1]
1 + ϕ[k −1]TP[k −1]ϕ[k −1]
.
(3.20)
Estas ecuaciones constituyen el núcleo del algoritmo RLS, el cual ajusta los parámetros del modelo en cada instante de muestreo utilizando la nueva información de entrada y salida. La ganancia adaptativa implícita en (3.19) actúa ponderando la magnitud del error y la información contenida en el vector regresor, mientras que la matriz P[k] regula la confianza en los parámetros actuales y controla la convergencia del proceso. En la práctica, el algoritmo puede incluir un factor de olvido λ ∈(0, 1] para dar mayor peso a los datos recientes y permitir una mejor adaptación frente a cambios en la dinámica del sistema. Este mecanismo es especialmente relevante en el cuadricóptero, donde las condiciones aerodinámicas y los efectos de acoplamiento pueden variar con el tiempo.
El esquema RLS implementado en este trabajo se ejecuta en paralelo con el controlador, actualizando los parámetros en tiempo real sin interferir con la estabilidad del lazo de control. De esta forma, el modelo identificado mantiene una correspondencia continua con la planta física, proporcionando al sistema adaptativo la información necesaria para el ajuste dinámico de las ganancias, tal como se desarrolla en el Capítulo 5. Esquema procedimental (RLS–ARX MIMO) Con el objetivo de mantener un modelo actualizado de la planta durante la operación, se implementa un esquema en línea basado en Mínimos Cuadrados Recursivos (RLS) aplicado a la estructura MIMO ARX. El flujo descansa en la formulación de regresión de la Subsección 3.1 ((3.9)–(3.10)) y en la ley de corrección paramétrica de (3.17)–(3.18), cuyo núcleo recursivo está dado por (3.19)–(3.20). A continuación se describen, de forma operativa y enlazada, la señal de excitación, la formación del regresor en tiempo real, el diagrama de bloques empleado y el pseudocódigo que ejecuta la actualización recursiva. Señales de excitación y preprocesamiento.
Para garantizar persistencia de excitación sin comprometer la seguridad de la planta, se emplea una excitación basada en señales binarias pseudoaleatorias (PRBS, por sus siglas en inglés Pseudo-Random Binary Sequence) en una configuración de múltiples entradas y múltiples salidas (MIMO), independiente por canal de entrada (q = 4). Se trabaja con un período de muestreo uniforme Ts = 0,01 s y una duración total de 50 s (N = 5000 muestras), utilizando una permanencia de 30 muestras por nivel y amplitudes por canal ampU = h
0,5, 0,05, 0,05, 0,05
i⊤.
CAPÍTULO 3. IDENTIFICACIÓN DE SISTEMAS MULTIVARIABLES
0
5
10
15
20
25
30
35
40
45
50
-0.5
0
0.5
u1 PRBS Excitation (u1 and u2)
0
5
10
15
20
25
30
35
40
45
50
Time [s]
-0.05
0
0.05
u2 (a) Señales de excitación PRBS para los canales u1 y u2.
0
5
10
15
20
25
30
35
40
45
50
-0.05
0
0.05
u3 PRBS Excitation (u3 and u4)
0
5
10
15
20
25
30
35
40
45
50
Time [s]
-0.05
0
0.05
u4 (b) Señales de excitación PRBS para los canales u3 y u4.
Figura 3.1. Señales de excitación PRBS en configuración MIMO por canal de entrada. La señal U ∈RN×4 se empaqueta como timeseries para su uso en el bloque From Workspace (variable U_ts) de MATLAB. Además, con el fin de emular condiciones realistas, las salidas medidas (p = 6) incorporan ruido Gaussiano aditivo de baja varianza, ymeas[k] = y[k] + v[k], con v[k] ∼N(0, σ2
y) y σy = 10−3. Finalmente, se considera un
período de inicialización warmup = m´ax(na, nb) + nk que asegura la disponibilidad de retardos para el primer vector regresor (3.10).
Esquema en Simulink y flujo de datos.
Sobre esa base, la Figura 3.2 sintetiza el flujo de identificación en línea. La señal U_ts alimenta, por un lado, la planta discreta cuadricóptero 6DOF, que entrega ymeas (con ruido), y, por otro, el bloque MakeDelaysMIMO, encargado de mantener memorias intermedias (buffers) circulares de y y u y de emitir los retardos Y_del y U_del. Con estos retardos, RLS_ARX_MIMO forma el vector ϕ[k −1] (3.10), calcula la predicción a un paso ˆy[k] = ˆθ[k −1]⊤ϕ[k −1] y actualiza ˆθ[k] y P[k] según (3.19)–(3.20). Las salidas yhat y theta_out se registran junto con ymeas. De forma opcional, la señal de reinicio reset reestablece ˆθ y P a ˆθ0 y P0 (p. ej., P0 = 106I), lo que permite reiniciar la estimación cuando es pertinente. Formación del regresor en tiempo real (MakeDelaysMIMO).
En línea con lo anterior, el bloque MakeDelaysMIMO organiza, para cada k, el vector ϕ[k −1] de acuerdo con (3.10) y la convención de signos de (3.11). El procedimiento se resume en el Algoritmo 2, cuyo objetivo es garantizar que la información más reciente de y y u alimente coherentemente el predictor.
Quadcopter 6DOF u y_meas Y_del U_del reset yhat Theta_out ready ++ y_meas u reset Y_del U_del ready Figura 3.2. Diagrama de bloques del esquema de identificación RLS–ARX MIMO en Simulink.
Algorithm 2 Formación de ϕ[k −1] en línea (bloque MakeDelaysMIMO) Require: na, nb, nk(= 1); secuencias ymeas[k], u[k]; k0 = m´ax(na, nb) + nk 1: Actualizar buffers circulares Y y U con ymeas[k] y u[k]. 2: if k ≥k0 then 3:
ϕ[k−1] ← h −y⊤[k−1], . . . , −y⊤[k−na], u⊤[k−nk], . . . , u⊤[k−nk−nb+1] i⊤ 4:
ready ←true 5: else 6:
ready ←false 7: end if 8: Emitir ϕ[k −1] y la bandera ready.
Actualización recursiva (RLS_ARX_MIMO).
En paralelo, el bloque RLS_ARX_MIMO realiza la actualización paramétrica con factor de olvido λ ∈(0, 1] y reinicio opcional, siguiendo la implementación vectorial de (3.19)–(3.20). El Algoritmo 3 recoge las operaciones principales y deja explícita la interacción entre la ganancia adaptativa, el error de estimación y la evolución de la covarianza, elementos que, en conjunto, sostienen la capacidad de adaptación del modelo durante la operación.
3.4.
Representación de estado extendida Con la estimación offline y online ya establecida, el siguiente paso consiste en dotar al modelo de una representación en espacio de estados que sea directamente utilizable para síntesis de control. Partimos de la forma de regresión introducida en el Capítulo, y[k] = Θ⊤ϕ[k −1],
(3.21)
CAPÍTULO 3. IDENTIFICACIÓN DE SISTEMAS MULTIVARIABLES
Algorithm 3 RLS–ARX MIMO con factor de olvido λ Require: λ ∈(0, 1], ˆθ0, P0 = αI (α ≫1); señales ymeas[k], ϕ[k −1]; reset 1: if reset then 2:
ˆθ ←ˆθ0, P ←P0 3: end if 4: if ready then 5:
ˆy[k] ←ˆθ⊤ϕ[k −1] 6:
e[k] ←ymeas[k] −ˆy[k] ▷Eq. (3.18) 7:
K[k] ← P ϕ[k −1] λ + ϕ[k −1]⊤P ϕ[k −1] 8:
ˆθ ←ˆθ + K[k] e[k] ▷Eq. (3.19) 9:
P ←1 λ P −K[k]ϕ[k −1]⊤P ▷Eq. (3.20) 10:
yhat ←ˆy[k], theta_out ←ˆθ 11: else 12:
Mantener ˆθ, P y yhat previos.
13: end if 14: Salvaguardas: simetrizar P, acotar cond(P), forzar 0 < λ ≤1. donde ϕ[·] concatena las salidas y entradas pasadas (hasta órdenes n, m) y Θ contiene los bloques −A1, . . . , −An, B1, . . . , Bm estimados. Definiendo el estado extendido como xa[k] := ϕ[k −1], se obtiene una realización discreta con estructura de registro de desplazamiento: xa[k + 1] = Fa[k] xa[k] + G u[k], y[k] = Ca[k] xa[k],
(3.22)
donde Ca[k] = Θ[k]⊤y las matrices Fa[k] y G adoptan la forma por bloques Fa[k] =
0
· · ·
0
0
0
· · ·
0
0
I
· · ·
0
0
0
· · ·
0
0
...
...
...
...
...
...
...
0
· · ·
I
0
0
· · ·
0
0
B1[k]
· · ·
Bm−1[k] Bm[k] −A1[k]
· · ·
−An−1[k] −An[k]
0
· · ·
0
0
I
0
· · ·
0
0
· · ·
0
0
0
...
...
...
0
· · ·
0
0
0
· · ·
I
0
,
G = I
0
...
0
.
(3.23)
La parte superior e inferior de Fa[k] implementa el desplazamiento de las muestras (bloques I y ceros), mientras que la fila media introduce el efecto de los coeficientes estimados
{Ai[k], Bi[k]}. Obsérvese que, cuando la identificación es online (Sección 3.3), Fa[k] y Ca[k] se actualizan con los parámetros corrientes, lo que produce un modelo dependiente del tiempo coherente con la dinámica observada.
Para su uso en control, esta realización se somete a chequeos estructurales ligeros: controlabilidad de (Fa[k], G) y observabilidad de (Fa[k], Ca[k]) mediante pruebas PBH en tiempo discreto, cuyos detalles se desarrollan en el Apéndice B (Sección B.3). En este trabajo, dichos chequeos actúan como condición de calidad antes de emplear el modelo en el rediseño de ganancias (Capítulo 5); si fallan, se bloquean las actualizaciones y se mantiene el modo nominal seguro.
3.4.1.
Control en espacio de estados discreto Considerando el modelo discreto x[k + 1] = F x[k] + G u[k],
(3.24)
y[k] = C x[k],
(3.25)
el diseño en espacio de estados requiere propiedades mínimas sobre el par (F, G) y el par (F, C). En particular, si el sistema es controlable, es posible asignar polos cerrados o resolver reguladores óptimos en tiempo discreto. La controlabilidad puede verificarse, por ejemplo, mediante la matriz W = h G FG · · · F n−1G i
,
rank(W) = n,
(3.26)
o de forma equivalente con el criterio PBH. En este trabajo, la verificación se aplica al modelo nominal obtenido offline y, cuando procede, al modelo extendido actualizado online; el diseño de control (pesos Q, R, acción integral) se desarrolla en detalle en el Capítulo 4, mientras que el rediseño adaptativo basado en la DARE discreta se presenta en el Capítulo 5.
En síntesis, la conversión ARX→espacio de estados establece una interfaz directa entre la identificación (que provee {Ai, Bi} y, por tanto, Fa, Ca) y el diseño de control discreto. Esta construcción mantiene la coherencia temporal del modelo con los datos y habilita, cuando las condiciones estructurales son satisfactorias, tanto el diseño nominal como las actualizaciones adaptativas de ganancias en operación. Como contraste con el enfoque recursivo adoptado, en el Apéndice A se presenta un estudio adicional basado en identificación por subespacios (MOESP/N4SID) sobre los mismos datos experimentales, lo que permite verificar el orden y la estructura del modelo discreto seleccionado.
CAPÍTULO 3. IDENTIFICACIÓN DE SISTEMAS MULTIVARIABLES
Con el modelo identificado en forma ARX, su actualización en línea mediante RLS y la realización discreta en espacio de estados (Fa, G, Ca), se dispone de una representación lista para síntesis de control. El Capítulo 4 desarrolla el diseño nominal en tiempo discreto (LQR + acción integral) apoyado en estas estructuras.
Capítulo 4 Diseño y Análisis del Controlador El diseño en tiempo discreto se apoya en una representación de la planta que admite realimentación de estados y tratamiento explícito de acoplamientos multivariables. El objetivo es garantizar la regulación y el seguimiento en régimen próximo al vuelo estacionario bajo restricciones físicas y computacionales, preservando la estabilidad interna y la coherencia con el ciclo de control; por ello, se integran desde el inicio límites de actuación, requisitos de referencia y latencias de ejecución.
Antes de sintetizar ganancias, se verifican condiciones discretas de controlabilidad y observabilidad como filtro operativo: cuando se satisfacen, habilitan la síntesis; en caso contrario, orientan ajustes mínimos de la realización o de la selección de salidas. Sobre esta base, las leyes de control se plantean mediante penalizaciones cuadráticas que equilibran esfuerzo y calidad de respuesta, incorporando acción integral para suprimir error estacionario y salvaguardas frente a saturación. El resultado es un controlador de ganancias fijas (controlador nominal de referencia) que actúa como ley estabilizante nominal y punto de retorno seguro, sobre el cual se articula el esquema adaptativo del Capítulo 5.
4.1.
Objetivos de control y restricciones El diseño se plantea en un régimen próximo al vuelo estacionario, donde la aproximación lineal local y los acoplamientos multivariables son dominantes. Se adoptan como salidas reguladas la posición (x, y, z) y el ángulo de guiñada ψ, mientras que el vector de control agrupa el empuje total y los momentos de actitud.
En este contexto, se busca una respuesta acotada y estable frente a referencias constantes o suavemente variables, supresión del error estacionario en las salidas seleccionadas y moderación del esfuerzo de control. En consecuencia, los pesos del regulador penalizan
CAPÍTULO 4. DISEÑO Y ANÁLISIS DEL CONTROLADOR
estados y entradas respetando la jerarquía de lazos (actitud más rápida que posición) y las restricciones de saturación; cuando procede se incorpora acción integral y estimación de estados en coherencia con la discretización y el nivel de ruido de medición. Aunque el diseño nominal asume disponibilidad de todos los estados, en el Apéndice D se analizan observadores tipo Luenberger, de orden completo y reducido, que ilustran cómo reconstruir estados a partir de un conjunto limitado de salidas. Tabla 4.1. Parámetros de diseño y restricciones operativas. Elemento Especificación Salidas reguladas / seguidas y = [ x y z ψ ]⊤ Entradas de control u = [ T τx τy τz ]⊤ Discretización Período de muestreo Ts = 0,01 s (100 Hz) Clases de referencia r[k]: señales constantes o multisinusoidales suaves Acción integral (LQI) Integradores en {x, y, z, ψ}: ξ[k+1] = ξ[k] + Ts (r[k] −y[k]) Pesos del regulador Q = Qx
0
0
Qi , R =
0,6
0
0
0
0
0,6
0
0
0
0
0,6
0
0
0
0
0,8
Qx = diag(3, 3, 6, 0,4, 0,4, 0,4, 0,2, 0,2, 0,6, 0,1, 0,1, 0,1) Qi =
40
0
0
0
0
40
0
0
0
0
60
0
0
0
0
50
Saturación de actuadores um´ın ≤u ≤um´ax, um´ax = [ 10, 0,6, 0,6, 0,4 ]⊤ um´ın = −um´ax Límite de tasa de cambio del mando (slew-rate)
∆u ∆t
∞≤˙um´ax Anti–windup Integración condicionada ante saturación Ruido de medición σ = [10−3, 10−3, 10−3, 5·10−4, 5·10−4, 5·10−4] Ry = diag(σ2) En conjunto, las especificaciones de la Tabla 4.1 fijan el entorno de diseño: el período de muestreo y la clase de referencias acotan la banda de interés; la estructura de Q y R prioriza la regulación de posición con esfuerzo moderado; la acción integral suprime el error estacionario en las salidas seleccionadas, y las restricciones de magnitud y tasa de cambio del mando, junto con Ry, condicionan la integración anti–windup y una estimación coherente con el nivel de ruido. Además del ajuste manual de Q y R, el Apéndice G
explora una sintonía automática mediante algoritmos genéticos sobre el mismo modelo y banco de pruebas; asimismo, una vez verificada la consistencia del modelo discreto mediante los chequeos estructurales de controlabilidad, detectabilidad y localización espectral documentados en el Apéndice B (Sección B.3) y resumidos en la Tabla B.3, se procede al diseño nominal en tiempo discreto.
4.2.
Diseño nominal en tiempo discreto El diseño se formula sobre una realización discreta (F, G, C) que satisface los chequeos estructurales previos y se ajusta a los parámetros y restricciones de la Tabla 4.1. Se adopta regulación lineal cuadrática y seguimiento mediante acción integral, en coherencia con Ts = 0,01 s, las salidas reguladas [x, y, z, ψ], las cotas de los actuadores y las ponderaciones Q y R allí fijadas.
4.2.1.
Regulador lineal cuadrático discreto (DLQR) Se considera la ley de realimentación u[k] = −K x[k], donde K se obtiene a partir de la DARE asociada al funcional cuadrático J = ∞ X k=0 x[k]⊤Q x[k] + u[k]⊤R u[k]
,
Q ⪰0, R ≻0.
Bajo las condiciones estándar de estabilizabilidad y detectabilidad, la ecuación de Riccati discreta P = F ⊤ P −P G R + G⊤PG −1G⊤P F + Q posee una solución simétrica definida positiva P, y la ganancia óptima viene dada por K = R + G⊤PG −1G⊤PF.
El escalado previo de x y u se mantiene coherente con las saturaciones de la Tabla 4.1 para evitar desbalanceos en Q y R y favorecer el acondicionamiento numérico. Como referencia, el Apéndice B recoge diseños alternativos por asignación espectral y modal de polos, útiles para comparación frente al enfoque óptimo LQR empleado.
CAPÍTULO 4. DISEÑO Y ANÁLISIS DEL CONTROLADOR
4.2.2.
Seguimiento con acción integral Para suprimir el error estacionario en [x, y, z, ψ], se introduce el estado integral del error ξ[k + 1] = ξ[k] + Ts r[k] −y[k]
,
y[k] = Csx[k], donde Cs selecciona las salidas reguladas. A partir de ello se define el sistema aumentado xa[k + 1] = F
0
−TsCs I | {z } Fa xa[k] + G
0
| {z } Ga u[k], xa = x ξ , y la ley de control u[k] = − h Kx Ki i xa[k], donde h Kx Ki i proviene de LQR sobre (Fa, Ga) con Q = blkdiag(Qx, Qi) y R según la Tabla 4.1. La estructura por bloques preserva la jerarquía de tiempos actitud–posición y la coherencia con las cotas de los actuadores.
Las respuestas detalladas del regulador de estado asociado, junto con métricas de desempeño para distintos escenarios, se documentan en el Apéndice C. El Apéndice F compara de forma sistemática el seguimiento de referencia mediante ganancia directa Kr y mediante acción integral, con métricas de error para cada salida, lo que respalda la elección de la estructura LQI (LQR + acción integral) adoptada.
4.2.3.
Saturaciones y anti–windup Las entradas se acotan en magnitud, um´ın ≤u[k] ≤um´ax, um´ax = h
10
0,6
0,6
0,4
i⊤, um´ın = −um´ax, y se restringe su cambio entre muestras para evitar maniobras bruscas, ∥u[k] −u[k −1] ∥∞≤∆um´ax.
Para la acción integral se aplica un esquema de anti–windup por retrocálculo, ξ[k + 1] = ξ[k] + Ts e[k] + Ts Kaw usat[k] −ucmd[k]
,
donde e[k] = r[k] −y[k], ucmd = −Kxx[k] −Kiξ[k] y usat[k] es la señal tras las limitaciones anteriores. Como salvaguarda adicional, puede congelarse la integración mientras exista saturación; en coherencia con Ts y los márgenes de los actuadores, Kaw se elige moderado para evitar oscilaciones.
(a) (b) Figura 4.1. Esquema general del controlador LQI. a) representación compacta del lazo integral con el regulador LQR; b) realización detallada en tiempo discreto en espacio de estados.
La Figura 4.1 muestra la implementación del controlador LQI. El esquema compacto muestra el lazo integral externo que corrige el error de seguimiento y alimenta al regulador LQR que gobierna la planta, mientras que el diagrama detallado explicita la realización discreta en espacio de estados. Para mantener la claridad, las saturaciones y el lazo de anti–windup no se representan de forma explícita, aunque se encuentran activas en la implementación.
4.3.
Control nominal ante incertidumbre y enlace adaptativo El modelo identificado constituye una base suficiente para el diseño, aunque sujeto a variaciones acotadas y ruido de medición. En consecuencia, la sintonía del controlador es conservadora: las ponderaciones de la Tabla 4.1 moderan el esfuerzo y refuerzan el amortiguamiento, la ubicación de los polos cerrados se mantiene ampliamente dentro del círculo unitario y se limita la variación de las señales de control entre muestras sucesivas. De este modo, pequeñas desviaciones de los coeficientes no comprometen la estabilidad ni la regularidad de la respuesta dentro del dominio de validez, y se preservan chequeos
CAPÍTULO 4. DISEÑO Y ANÁLISIS DEL CONTROLADOR
estructurales en instantáneas representativas para verificar que las condiciones mínimas continúan satisfechas.
Sobre esta base, el controlador de ganancias fijas garantiza operación segura y actúa como soporte del esquema adaptativo. El sistema expone indicadores de calidad de modelo y de disponibilidad de actuadores; cuando son favorables, el ajuste de ganancias se habilita de forma gradual y con tiempo mínimo de permanencia, mientras que, si se degradan, el retorno al controlador base es inmediato. La convivencia entre diseño nominal y actualización en línea se realiza así con transiciones suaves y bajo protección frente a saturaciones y cambios bruscos de mando, estableciendo un marco robusto sobre el cual se desarrolla el esquema adaptativo del capítulo siguiente.
Capítulo 5 Control Adaptativo Óptimo: STR Indirecto En este capítulo se presenta un marco de control adaptativo de tipo regulador autoajustable indirecto (STR, por sus siglas en inglés Self-Tuning Regulator) para el cuadricóptero, en el que la identificación recursiva en línea mediante mínimos cuadrados recursivos (RLS) se combina con el rediseño periódico de las ganancias de un regulador óptimo discreto con acción integral. Este rediseño se acompaña con una ganancia de referencia Kr calculada sobre el modelo identificado, de modo que la ley adaptativa conserve la estructura LQI y, al mismo tiempo, incorpore una compensación anticipativa de referencia. El objetivo es que el sistema se ajuste de forma automática ante variaciones paramétricas, perturbaciones exógenas y efectos no modelados, preservando la estabilidad interna y garantizando un desempeño robusto dentro del dominio operativo de interés. Como complemento, el Apéndice H presenta una extensión LQG que integra el LQR con un estimador de Kalman en tiempo discreto y resultados de seguimiento en presencia de ruido significativo en los sensores.
La propuesta se articula sobre tres pilares complementarios. En primer lugar, se asume un modelo nominal en tiempo discreto, derivado del análisis del Capítulo 2.6 y refinado, cuando es necesario, mediante los procedimientos de identificación del Capítulo 3. En segundo lugar, se implementa un estimador paramétrico en línea que actualiza un modelo entrada–salida (por ejemplo, ARX-MIMO) a partir de datos de operación, garantizando excitación persistente y control de la variabilidad mediante regularización y olvido exponencial. En tercer lugar, se incorpora un mecanismo de autoajuste que, a intervalos programados o bajo criterios de activación basados en métricas de calidad del modelo, resuelve la ecuación de Riccati discreta correspondiente y reconfigura adaptativamente las
CAPÍTULO 5. CONTROL ADAPTATIVO ÓPTIMO: STR INDIRECTO
ganancias del regulador óptimo, respetando restricciones prácticas como saturaciones de actuadores, límites de esfuerzo y tiempos mínimos de permanencia para evitar conmutaciones rápidas.
Como continuación del Capítulo 4, donde se estableció el diseño nominal de controladores óptimos en tiempo discreto bajo supuestos de modelo fijo, este capítulo aborda un escenario operativo en el que los parámetros de la planta pueden variar y las incertidumbres adquieren un papel central. En este contexto, el énfasis se sitúa en integrar de manera coherente la identificación y el control óptimo dentro de un ciclo de adaptación: a partir de datos en línea se actualiza el modelo mediante técnicas recursivas y, con base en estas actualizaciones, se rediseñan periódicamente las ganancias del regulador, cerrando así un bucle Modelado →Identificación →Control Adaptativo que preserva la estabilidad y mejora el desempeño frente a perturbaciones y deriva paramétrica (Figura 5.1). Esta perspectiva consolida el puente metodológico entre los capítulos previos de modelado e identificación y la síntesis de control, y sienta las bases para la evaluación experimental mediante simulación.
Desde el punto de vista teórico, el STR indirecto explora el principio de separación en su versión discreta: la ley de control por realimentación de estados (LQR/LQI) y el estimador de estado (Luenberger/Kalman) pueden diseñarse de forma desacoplada siempre que el par identificado cumpla condiciones de estabilizabilidad y detectabilidad. Para escenarios con ruido de medición explícito, el Apéndice E desarrolla la formulación del estimador óptimo en el sentido de Kalman y detalla el diseño de filtros LQE coherentes con el marco de control considerado. Sobre esta base se formalizan garantías de estabilidad local bajo actualización de parámetros, se analiza la interacción entre errores de identificación y síntesis de ganancias y se discute la sensibilidad frente a ruido de medición y perturbaciones no estacionarias, incorporando además consideraciones de implementación embebida relacionadas con el costo computacional, las limitaciones de hardware y la robustez numérica en el cálculo en tiempo real.
En síntesis, el contenido que sigue establece el andamiaje matemático del STR indirecto, describe su realización algorítmica y propone criterios operativos para su uso seguro en lazo cerrado. Ello prepara el terreno para su integración con los módulos de percepción y planificación, así como para su evaluación en simulación en los capítulos posteriores.
Modelado (Cap. 2) Identificación (Cap. 3) Control Óptimo (Cap. 4) Control Adaptativo STR Indirecto (Cap. 5) {A, B, C, D} nominal bΘ, \ {A, B, C, D} Kx, Ki iniciales datos en línea (u, y) modelo actualizado / Kx, Ki, Kr reajustados Figura 5.1. Esquema conceptual del flujo Modelado →Identificación →Control Óptimo →Control Adaptativo (STR indirecto). Las flechas sólidas representan el flujo nominal entre capítulos; las flechas discontinuas ilustran el bucle de adaptación con datos en línea y rediseño periódico de ganancias.
5.1.
Marco General del STR Indirecto El controlador adaptativo se basa en un principio sencillo: mantener la estructura nominal fijada en el Capítulo 4 y reajustar sus ganancias cuando el modelo actualizado aporta información suficiente y consistente. Para ello, el esquema adopta un enfoque indirecto en el que se actualiza en línea un modelo discreto de la planta y, en instantes intermitentes, se recalculan las ganancias de un LQI junto con una ganancia de referencia Kr, empleando las ponderaciones y restricciones ya definidas. De este modo, la adaptación no altera el ciclo de control, respeta el periodo de muestreo, conserva las salidas reguladas [x, y, z, ψ] y mantiene explícitos los límites de magnitud y de cambio entre muestras. La operación inicia con un LQI nominal físico de 12 estados, usado como modo seguro y como respaldo permanente durante toda la simulación.
Este planteamiento contrasta con el enfoque directo, donde las ganancias se ajustan de manera inmediata a partir del error. En el enfoque indirecto se separan identificación y control: el primer bloque estima el comportamiento de la planta y el segundo resuelve un problema lineal cuadrático discreto sobre el modelo actualizado. Esta organización se ilustra en la Figura 5.2, donde el regulador autoajustable envuelve el lazo de control nominal e integra los módulos de identificación y rediseño del controlador. Así se preserva la compatibilidad con el marco óptimo en tiempo discreto, se facilita la previsualización del mando frente a saturaciones y se mantiene la coherencia con la jerarquía actitud–posición ya definida; cuando la evidencia del identificador es insuficiente o los actuadores se hallan cerca de sus límites, las ganancias no se modifican y el controlador nominal continúa
CAPÍTULO 5. CONTROL ADAPTATIVO ÓPTIMO: STR INDIRECTO
gobernando.
El flujo operativo es regular. En cada muestra se adquieren las señales, se actualizan los parámetros del modelo mediante un estimador recursivo con factor de olvido y se registran indicadores de calidad suavizados. Periódicamente, el modelo se convierte a una realización discreta en espacio de estados alineada con las salidas de interés. Con esta información, un supervisor decide si procede un nuevo cálculo de [ Kx Ki ] y de la ganancia de referencia Kr: sólo cuando los indicadores resultan favorables y la demanda de actuadores se mantiene dentro de márgenes se resuelve la Riccati discreta, se calcula la compensación anticipativa de referencia y se verifica, por adelantado, que el mando previsto no vulnera ni la magnitud ni la variación admisible. En tal caso, las ganancias se aplican con mezcla gradual para evitar transitorios bruscos; en caso contrario, se difiere el reajuste. En todo momento permanecen activas las salvaguardas de saturación y el anti-windup discreto de la Subsección 4.2.3, de modo que la adaptación convive con las mismas protecciones que rigen el diseño nominal.
En conjunto, el STR indirecto actúa como un puente entre el controlador de ganancias fijas y la actualización en línea: cuando la información lo permite, ajusta; cuando no, sostiene la operación con el diseño base. Esta lógica, apoyada en métricas estables y en decisiones con histéresis y tiempo mínimo de permanencia, proporciona transiciones suaves y una integración natural con el marco de control óptimo discreto. Figura 5.2. Arquitectura general de un regulador autoajustable (STR) indirecto. El módulo de identificación de sistema estima los parámetros del proceso y el bloque de diseño actualiza en línea los parámetros del controlador que actúa sobre la planta.
5.2.
Sistema de Identificación El lazo adaptativo descansa en un estimador en línea que, muestra a muestra, actualiza un modelo discreto alineado con el régimen de operación (vuelo estacionario). La relación entre entrada y salida se expresa como regresión (Subsección 3.1) y los parámetros se estiman mediante mínimos cuadrados recursivos con factor de olvido (Sección 3.3.1). La actualización del estimador opera al periodo de muestreo Ts = 0,01 s, en coherencia con la discretización adoptada en el Capítulo 4. En la implementación, el identificador ARX-MIMO se concentra en los canales regulados [x, y, z, ψ], coherentes con las salidas evaluadas y con la referencia de seguimiento. En cada muestra se arma el vector de regresores con pasados de estas salidas y de las entradas de control, se actualizan los parámetros y, cuando es necesario, se aplica una proyección suave para acotar cambios improbables. El identificador procesa las mediciones, calcula el error de predicción y registra indicadores de calidad (ajuste suavizado, condición numérica y correlaciones con la entrada) que otorgan mayor peso a la información reciente y que el supervisor empleará más adelante para habilitar o posponer cualquier reajuste del controlador. Cuando la excitación es insuficiente, puede añadirse una señal de prueba de baja amplitud, compatible con las restricciones de la Tabla 4.1, de modo que se sostenga el aprendizaje sin perturbar la tarea principal. Con ello, el modelo se mantiene alineado con la operación y acumula métricas de confianza que el supervisor usará para tomar decisiones de reajuste. Con los parámetros actualizados, en instantes de decisión prefijados se convierte el modelo ARX-MIMO a una realización en espacio de estados tipo shift-register. Esta representación no corresponde al estado físico completo del cuadricóptero, sino a un estado interno asociado a la memoria del modelo entrada–salida, y permite disponer de las matrices (Fa[k], Ga, Ca[k]) requeridas para el rediseño del controlador. Este paso, presentado en el Capítulo 3, deja preparado el modelo para recalcular, cuando corresponda, las ganancias del LQI y la ganancia de referencia Kr con las ponderaciones de la Tabla 4.1. En consecuencia, el identificador no solo actualiza parámetros, sino que entrega al supervisor una realización de diseño junto con sus métricas de confianza, de modo que se pueda decidir entre retocar las ganancias o mantener el controlador nominal.
CAPÍTULO 5. CONTROL ADAPTATIVO ÓPTIMO: STR INDIRECTO
5.3.
Reajuste discreto del controlador sobre el modelo actualizado El reajuste se ejecuta en instantes de decisión espaciados respecto del periodo de muestreo Ts, con el fin de acotar la carga computacional. En cada uno de estos instantes, a partir de las instantáneas (Fa[k], Ga, Ca[k]) y de las ponderaciones Q, R de la Tabla 4.1, se resuelve la DARE y se obtiene un candidato para las ganancias del LQI, h Ktry x Ktry i i
.
Sobre el mismo modelo identificado se calcula además una ganancia de referencia Ktry r , orientada a mejorar el seguimiento de referencia sin modificar la estructura integral del controlador.
Antes de adoptar el candidato, se previsualiza el mando asociado: uprev = −Ktry x xsr[k] −Ktry i ξ[k] + Ktry r r[k], donde xsr representa el estado asociado a la realización del modelo ARX y ξ[k] el estado integral del LQI. Para evitar saltos durante la conmutación, el estado integral se compatibiliza con el mando aplicado antes de aceptar el nuevo conjunto de ganancias. Posteriormente, se verifica que el mando previsto respete las cotas de magnitud y los límites de variación entre muestras establecidos en la Subsección 4.2.3; además, se exige que las métricas del identificador se mantengan dentro de umbrales de calidad durante una ventana mínima.
Cuando las comprobaciones son favorables, la transición se aplica por mezcla gradual, Kx ←(1 −αx)Kx + αxKtry x , Ki ←(1 −αi)Ki + αiKtry i
,
Kr ←(1 −αr)Kr + αrKtry r , lo que suaviza el cambio de realimentación, acción integral y compensación anticipativa, respetando la disponibilidad de actuadores. En caso contrario, la actualización se difiere hasta el siguiente instante de decisión. El proceso se registra con marcas de tiempo y referencias para facilitar el seguimiento del historial de ajustes.
5.4.
Gestión adaptativa En el esquema propuesto, el supervisor coordina el intercambio entre el identificador en línea, el modo seguro nominal y el controlador adaptativo. A intervalos regulares recibe instantáneas del modelo y sus indicadores de calidad, los contrasta con la disponibilidad
Algorithm 4 Reajuste discreto y mezcla suave de ganancias (STR indirecto) Require: Fa[k], Ga, Ca[k], Q, R, xsr(k), ξ(k), r(k); ganancias actuales Kx, Ki, Kr; límites um´ın, um´ax, ∆um´ax; banderas calidad_ok; instantes k, kult; parámetros αx, αi, αr ∈(0, 1), ventana mínima Nm´ın 1: if calidad_ok ∧(k −kult) ≥Nm´ın then 2:
P ←Riccati(Fa[k], Ga, Q, R) 3:
[Ktry x , Ktry i ] ←DLQI(Fa[k], Ga, Ca[k], Q, R) 4:
Ktry r ←GananciaReferencia(Fa[k], Ga, Ca[k], Ktry x ) 5:
uprev ←−Ktry x xsr(k) −Ktry i ξ(k) + Ktry r r(k) ▷Previsualización 6:
if DentroDeCotas(uprev, um´ın, um´ax) ∧VariacionAceptable(uprev, u[k−1], ∆um´ax) then 7:
Kx ←(1 −αx) Kx + αx Ktry x 8:
Ki ←(1 −αi) Ki + αi Ktry i 9:
Kr ←(1 −αr) Kr + αr Ktry r 10:
kult ←k 11:
end if 12: end if de actuadores y con las restricciones de la Tabla 4.1 y decide si activar el ajuste descrito en la Sección 5.3 o mantener (o restituir) el modo nominal.
En coherencia con el enfoque adoptado, la política de supervisión es deliberadamente conservadora: se exige una ventana en la que las métricas permanezcan dentro de umbrales y, además, se respeten los límites de magnitud y de variación entre muestras fijados en la Subsección 4.2.3. Se procesan indicadores compactos, tales como el ajuste suavizado del identificador, la condición numérica, medidas de blancura y la correlación entrada–error, así como la fracción de tiempo en saturación o próxima a los límites. Si estos indicadores se sostienen durante Ton, se habilita el modo adaptativo; si alguno se degrada de forma persistente durante Toff, se retorna al nominal. Adicionalmente, se impone un tiempo mínimo de permanencia Tdwell y un breve periodo de recuperación tras cada retorno. Asimismo, los reajustes se bloquean temporalmente ante saturación reciente o durante ventanas de protección alrededor de cambios de referencia, evitando rediseños en instantes dominados por transitorios fuertes.
La conmutación se ejecuta mediante mezcla suave de ganancias, empleando factores diferenciados para la realimentación, la acción integral y la ganancia de referencia, denotados por (αx, αi, αr), respectivamente. Esta transición atenúa saltos y preserva márgenes de los actuadores. Si se detectan condiciones de riesgo, la autorización se revoca de in-
CAPÍTULO 5. CONTROL ADAPTATIVO ÓPTIMO: STR INDIRECTO
Algorithm 5 Máquina de estados del supervisor Require: estado inicial ∈{NOMINAL, ADAPTATIVO}; funciones OK(Ton), Falla(Toff), SaturacionProlongada(), Riesgo(); manejadores IniciarDwell(), IniciarCooldown(); prueba DwellCumplido() 1: if estado = ADAPTATIVO ∧DwellCumplido() ∧Falla(Toff) then 2:
estado ←NOMINAL 3:
IniciarCooldown() 4: end if 5: if estado = ADAPTATIVO ∧Riesgo() then 6:
estado ←NOMINAL 7:
IniciarCooldown() 8: end if 9: if estado =
ADAPTATIVO
∧ DwellCumplido() ∧ Falla(Toff) ∨ SaturacionProlongada() then 10:
estado ←NOMINAL 11:
IniciarCooldown() 12: end if mediato y se mantiene o restituye el modo nominal. En este contexto, Riesgo() agrupa condiciones como saturación persistente, error de seguimiento excesivo o actitud fuera de margen. Todas las decisiones se registran con marcas de tiempo e indicadores asociados. Durante el Reajuste y la conmutación permanecen activos los límites de magnitud y de variación de u, así como el esquema de anti-windup definido en el Capítulo 4; en consecuencia, cualquier actualización de [ Kx Ki ] y de la ganancia de referencia Kr respeta dichos mecanismos.
Con el esquema adaptativo ya delineado queda cerrado el ciclo entre modelo, datos y control. En continuidad, el capítulo siguiente aborda la navegación autónoma apoyada en información espacial del entorno, en la que el seguimiento de las referencias de trayectoria descansa en el controlador adaptativo desarrollado en este capítulo. De este modo, las decisiones de planificación de rutas y replanificación se integran de forma coherente con el identificador y el supervisor, manteniendo los marcos temporales y las restricciones establecidos en la fase de diseño.
Capítulo 6 Percepción visual para navegación autónoma La navegación autónoma en entornos tridimensionales exige que el cuadricóptero disponga de una representación interna del espacio circundante, capaz de reflejar obstáculos y regiones libres de forma compatible con los algoritmos de planificación y control. En este contexto, la percepción visual se concibe como el conjunto de procesos que transforman información geométrica del entorno en una representación voxelizada de ocupación, actualizada a partir de datos sensoriales e integrada con los módulos de planificación global, planificación local y seguimiento de trayectorias, sobre los que se apoya el esquema de control adaptativo desarrollado en el Capítulo 5 para el seguimiento de referencias. Para materializar esta representación, se integró la librería map_manager [65], que a partir de la información de odometría y de la cámara de profundidad construye el mapa de ocupación voxelizado empleado por los algoritmos de planificación y por el esquema de control de navegación. Esta representación espacial discretizada dota al vehículo de capacidad para anticipar conflictos y adaptar el movimiento, y se consolida como un componente estructural de la arquitectura de navegación autónoma que enlaza los datos de simulación con las decisiones de planificación, control y supervisión.
6.1.
Arquitectura general de percepción y navegación basada en mapas de ocupación La arquitectura de navegación basada en percepción visual se implementa sobre Robot Operating System (ROS) Melodic, ejecutado en una plataforma embebida NVIDIA Jetson Nano con Ubuntu 18.04, que aloja los nodos de percepción, planificación y seguimiento
CAPÍTULO 6. PERCEPCIÓN VISUAL PARA NAVEGACIÓN AUTÓNOMA
de trayectorias. La dinámica del cuadricóptero y el entorno tridimensional tipo laboratorio se simulan en Gazebo, mientras que RViz, la herramienta de visualización tridimensional de ROS, se emplea para representar el estado del vehículo, el mapa de ocupación y las trayectorias. El flujo de información se organiza alrededor de un mapa de ocupación tridimensional, que constituye la representación común sobre la que actúan el planificador global, el módulo de suavizado de trayectorias, el seguidor de trayectorias y el planificador local reactivo.
La estructura global de este flujo se resume en la Figura 6.1, donde se aprecia cómo la nube de puntos del entorno se transforma en un mapa voxelizado, es decir, en una discretización tridimensional del espacio basada en vóxeles (voxels), que alimenta al planificador global basado en árboles aleatorios de exploración rápida (RRT, por sus siglas en inglés Rapidly-exploring Random Tree), al suavizador mediante curvas B-spline (Basis spline) y al planificador local, mientras que el seguidor de trayectorias genera referencias de pose para el controlador a partir de la trayectoria suavizada.
Sistema de Percepción y Navegación Basado en Voxel Map mediciones Sensores a bordo
(RGB-D / IMU)
Nube de puntos / Medida de Profundidad mapa 3D global VoxelGrid3D (mapa de ocupación voxelizado) trayectoria bruta Global Planner
RRT-3D
trayectoria suave B-Spline Smoother (suavizado global) parámetros / setpoints PathFollower3D (seguimiento de trayectoria 3D) mapa 3D local / referencia global trayectoria local LocalPlanner3D (planificación local reactiva y replanning sobre el mapa) señales a motores Controlador Quad (UAV / dron) Figura 6.1. Arquitectura general del subsistema de percepción y navegación basada en mapas de ocupación, implementada en ROS Melodic sobre una NVIDIA Jetson Nano. El mapa voxelizado se comparte entre los módulos de planificación global, suavizado de trayectorias, seguimiento y planificación local.
La Tabla 6.1 recoge los nodos ROS que materializan esta arquitectura y su función dentro del flujo de navegación. En las secciones siguientes se describen con mayor detalle sus interfaces y el modo en que interactúan con el mapa de ocupación y con el controlador. Tabla 6.1. Nodos ROS principales del subsistema de percepción y navegación. Nodo Función principal rrt_3d_interactive Planificador global RRT con interfaz RViz.
bspline_smoother Suavizado de rutas RRT mediante B-splines.
path_follower_3d Seguimiento de trayectorias 3D y gestión de misión.
local_planner_3d Replanificación local reactiva y costura de trayectorias.
6.2.
Entorno de simulación y flujo de ejecución El sistema de navegación se valida en un entorno de simulación compuesto por Gazebo y RViz sobre ROS Melodic. Gazebo modela la dinámica del cuadricóptero y un escenario interior configurado como un recinto cerrado con múltiples obstáculos cúbicos, mientras que RViz se utiliza para visualizar el mapa de ocupación voxelizado, la posición estimada del vehículo y las trayectorias generadas por los planificadores. La Figura 6.2 ilustra esta configuración: a la izquierda se muestra el mapa voxel coloreado según la altura y una trayectoria planificada, y a la derecha la escena equivalente en Gazebo con la disposición de obstáculos en planta.
El flujo de ejecución parte del lanzamiento conjunto de Gazebo y RViz mediante archivos launch que cargan el mundo del laboratorio, el modelo del cuadricóptero y la configuración de visualización. A continuación se activa el flujo de mapeo tridimensional, que obtiene nubes de puntos del escenario, construye el mapa de ocupación y lo publica en un tópico dedicado. Sobre este mapa opera el planificador global basado en RRT, cuya salida se suaviza mediante B-splines y se entrega al seguidor de trayectorias, mientras que el planificador local supervisa en paralelo la trayectoria y el mapa para detectar posibles conflictos y disparar replanificaciones cuando es necesario.
La coordinación entre estos nodos se apoya en archivos de configuración en formato YAML que definen parámetros geométricos del entorno, resoluciones de mapa, márgenes de seguridad y características de los algoritmos de planificación y seguimiento, lo que permite ajustar el comportamiento del sistema sin modificar el código fuente y reproducir de forma consistente el mismo flujo de ejecución en distintas sesiones de simulación.
CAPÍTULO 6. PERCEPCIÓN VISUAL PARA NAVEGACIÓN AUTÓNOMA
Figura 6.2. Vista combinada del entorno de simulación en Gazebo y su representación voxelizada en RViz. A la derecha se observa la disposición física de los obstáculos, mientras que a la izquierda se muestra el mapa 3D discretizado utilizado por los módulos de planificación y navegación.
6.3.
Representación voxelizada del entorno
6.3.1.
Definición y construcción de la grilla de ocupación La representación geométrica que sustenta la planificación global y la replanificación local se organiza en torno a una malla de vóxeles VoxelGrid3D, que discretiza el espacio tridimensional en celdas cúbicas de resolución uniforme. Esta estructura se define a partir de un origen origin_x, origin_y, origin_z, de un número de celdas por eje size_x, size_y, size_z y de una resolución espacial map_resolution que fija la longitud de arista de cada voxel. Cada celda se asocia a un estado de ocupación que permite distinguir entre regiones libres, ocupadas y exteriores al dominio considerado, y se actualiza a partir del mapa tridimensional mantenido por la librería map_manager [65].
El paso de coordenadas continuas (x, y, z) a índices discretos (ix, iy, iz) se realiza mediante ix = x −origin_x map_resolution
,
iy = y −origin_y map_resolution
,
iz = z −origin_z map_resolution
,
(6.1)
de forma que cada punto del espacio se proyecta sobre una única celda de la grilla. La elección de map_resolution establece un compromiso entre detalle geométrico y coste computacional: valores pequeños capturan mejor la forma de los obstáculos, pero incrementan el tamaño de la grilla y el tiempo requerido para planificación y comprobación de
colisiones.
La información necesaria para poblar la grilla procede de nubes de puntos tridimensionales publicadas en ROS como mensajes sensor_msgs/PointCloud2. Un nodo de mapeo recorre cada punto (x, y, z) de la nube, aplica la transformación de la ec. (6.1) y, siempre que los índices resultantes se encuentren dentro de los límites de la grilla, marca la celda correspondiente como ocupada; los puntos que caen fuera de dichos límites se descartan para evitar accesos fuera de rango. Este proceso transforma una nube de puntos densa en una representación binaria o ternaria de ocupación que preserva la geometría esencial del entorno sin almacenar todas las muestras individuales. El mapa resultante se publica periódicamente en un tópico, /occupancy_map/voxel_map, que alimenta los módulos de planificación global y local.
Figura 6.3. Discretización de una nube de puntos tridimensional en una malla de vóxeles con resolución fija. Cada punto se proyecta sobre una celda de la grilla, que se marca como ocupada si pertenece al dominio del mapa.
6.3.2.
Mapas voxel precomputados Además del mapeo en línea, se contempla el uso de mapas de ocupación precomputados almacenados en archivos .npz. Estos archivos contienen arreglos de puntos tridimensionales que describen un entorno fijo. Durante la inicialización, el nodo de planificación global puede cargar el archivo, determinar una caja envolvente (bounding box) del conjunto de puntos o emplear límites proporcionados por el usuario, e instanciar la estructura VoxelGrid3D a partir de dichos límites.
CAPÍTULO 6. PERCEPCIÓN VISUAL PARA NAVEGACIÓN AUTÓNOMA
Una vez creada la grilla, los puntos almacenados se recorren y se proyectan a índices discretos mediante (6.1), marcando como ocupadas las celdas correspondientes. El resultado es un mapa voxel estático que se publica en el mismo tópico que el mapa generado en línea, de forma que los módulos de planificación y seguimiento operan de manera transparente respecto al origen de la información. Esta conmutación entre mapas generados y precomputados es especialmente útil para la depuración del sistema y para campañas de experimentos reproducibles, en las que se mantiene fija la geometría del entorno mientras se modifican parámetros de planificación o de control.
6.4.
Planificación global sobre mapas de ocupación tridimensionales
6.4.1.
Nodo interactivo de planificación RRT La planificación global de trayectorias se realiza mediante un algoritmo Rapidlyexploring Random Trees (RRT) formulado sobre el mapa de ocupación tridimensional descrito en la Sección 6.3.
Este algoritmo se implementa en el nodo rrt_3d_interactive_node.py, que actúa como interfaz entre la representación voxelizada del entorno y el operador en RViz. El nodo suscribe el tópico /move_base_simple/goal, de modo que dos selecciones consecutivas con la herramienta 2D Nav Goal definen la posición de partida y la posición objetivo de la misión en el plano de vuelo.
Las posiciones seleccionadas se proyectan sobre celdas libres del mapa mediante funciones de tipo project_to_nearest_free, garantizando que tanto el punto inicial como el objetivo se encuentren en regiones transitables. A partir de estos puntos y del mapa VoxelGrid3D, el nodo instancia un objeto RRTPlanner2p5D que explora el espacio libre en el plano horizontal a una altura de referencia z_ref. Cuando se encuentra un camino que conecta inicio y objetivo sin penetrar en vóxeles ocupados, se construye un mensaje nav_msgs/Path con la secuencia de poses (x, y, zref) y se publica en el tópico /dynamicNavigation/rrt_path, desde donde es utilizado por el módulo de suavizado.
La Figura 6.4 muestra un ejemplo de ruta generada por este nodo sobre el mapa de ocupación del laboratorio simulado. El camino resultante aparece como una polilínea roja que recorre las regiones libres entre los bloques que representan los obstáculos, evidenciando la capacidad del planificador para explotar la estructura del mapa voxelizado y
producir trayectorias globales libres de colisiones.
Figura 6.4. Ejemplo de ruta global generada por el planificador RRT sobre el mapa de ocupación voxelizado del laboratorio. La trayectoria poligonal en rojo conecta el punto inicial con el objetivo evitando las regiones ocupadas representadas por los bloques azules.
6.4.2.
Formulación geométrica del RRT El planificador opera en el plano horizontal sobre el espacio libre Xfree ⊂R2, definido como el conjunto de puntos cuya proyección tridimensional no atraviesa celdas ocupadas de la grilla VoxelGrid3D. El algoritmo construye incrementalmente un árbol enraizado en xstart y lo expande con muestras aleatorias xrand ∈Xfree. Para cada muestra se selecciona el nodo más cercano xnear (distancia euclídea) y se propone xnew = steer xnear, xrand, δ
,
(6.2)
donde steer(·) avanza desde xnear hacia xrand con un paso máximo δ > 0. El tramo [xnear, xnew] se incorpora solo si permanece dentro de Xfree, verificado mediante consultas de ocupación sobre el mapa voxelizado.
La construcción del árbol continúa hasta que la distancia entre el nodo más próximo al objetivo y la posición deseada es inferior a una tolerancia goal_reach_distance, o bien hasta que se alcanza un tiempo máximo de planificación timeout. En el primer caso se obtiene una ruta poligonal que conecta inicio y objetivo recorriendo los nodos del árbol; en el segundo, el intento de planificación se considera fallido.
CAPÍTULO 6. PERCEPCIÓN VISUAL PARA NAVEGACIÓN AUTÓNOMA
6.4.3.
Configuración mediante parámetros en YAML El comportamiento del planificador global se controla mediante parámetros almacenados en archivos YAML que se cargan en el servidor de parámetros de ROS durante el lanzamiento del sistema. Esta configuración incluye, entre otros, el tamaño de la caja de colisión collision_box que aproxima el volumen del cuadricóptero, los límites del área de planificación env_box, la distancia de expansión rrt_incremental_distance, la probabilidad de conexión directa al objetivo rrt_connect_goal_ratio, la tolerancia de llegada goal_reach_distance y el tiempo máximo de búsqueda timeout. La Tabla 6.2 recoge algunos de estos parámetros representativos junto con su interpretación y valores típicos en el entorno de laboratorio.
Esta separación entre lógica algorítmica y configuración numérica permite ajustar el compromiso entre calidad de las trayectorias, robustez frente a obstáculos y coste computacional sin modificar el código fuente. Variaciones en el tamaño de la caja de colisión o en la resolución de la grilla se reflejan en el margen de seguridad que el planificador mantiene respecto a paredes y obstáculos interiores, mientras que la elección de δ, timeout y rrt_connect_goal_ratio condiciona el tiempo de respuesta y la probabilidad de éxito de la planificación.
Tabla 6.2. Parámetros representativos del planificador global RRT. Todos los parámetros corresponden a la configuración del módulo RRT global.
Parámetro Significado Valor típico collision_box Caja de colisión del cuadricóptero [0.4, 0.4, 0.2] m rrt_incremental_distance Distancia máxima de expansión δ
1.0 m
rrt_connect_goal_ratio Probabilidad de conexión al objetivo
0.5
goal_reach_distance Tolerancia de llegada al objetivo
0.5 m
6.5.
Suavizado de trayectorias mediante B-splines Las trayectorias calculadas por el planificador global basado en RRT se obtienen como una sucesión de segmentos rectos que conectan los nodos del árbol en el espacio libre. Esta construcción poligonal garantiza la factibilidad geométrica del camino, pero no impone restricciones sobre la curvatura ni sobre la variación de la dirección entre segmentos consecutivos, dando lugar a trayectorias con quiebres pronunciados que pueden traducirse en demandas abruptas sobre los actuadores y en requerimientos de aceleración poco com-
patibles con un vuelo suave. En un vehículo aéreo multirrotor, estos cambios bruscos de rumbo o de velocidad de referencia incrementan el esfuerzo del controlador y dificultan el uso de la ruta como referencia continua para estrategias de control avanzadas o para módulos de planificación local que asumen cierto grado de regularidad geométrica. El suavizado mediante B-splines ofrece un compromiso entre fidelidad a la ruta global y mejora de las propiedades geométricas de la trayectoria. A partir de la secuencia de puntos generada por el planificador, se construye una curva suave que mantiene la estructura general del camino, pero redistribuye la curvatura y elimina esquinas pronunciadas. Esta curva suavizada se utiliza como referencia para el seguidor de trayectorias, reduciendo las variaciones abruptas de orientación y facilitando un comportamiento más regular del cuadricóptero.
6.5.1.
Nodo de suavizado y algoritmo de Chaikin El suavizado de la trayectoria global se implementa en un nodo dedicado de ROS, que suscribe el tópico /dynamicNavigation/rrt_path y publica la trayectoria procesada en /dynamicNavigation/bspline_trajectory. El nodo extrae la secuencia de puntos bidimensionales {pi} asociada a la ruta RRT y aplica iterativamente el algoritmo de Chaikin, utilizando como parámetros principales el número de iteraciones iterations y la distancia mínima entre puntos consecutivos min_segment_length, que controlan el grado de suavizado y la densidad final de la trayectoria. En cada iteración, para cada par consecutivo de puntos (pi, pi+1) se generan dos nuevos puntos qi y ri definidos como qi = 0,75 pi + 0,25 pi+1, ri = 0,25 pi + 0,75 pi+1,
(6.3)
mientras que los puntos extremos de la secuencia original se conservan. La repetición del proceso produce una curva que aproxima una B-spline cuadrática, suavizando progresivamente los quiebres presentes en la ruta inicial. Una vez completadas las iteraciones, se realiza un remuestreo espacial que elimina puntos demasiado próximos entre sí, de forma que la distancia entre puntos consecutivos sea superior al umbral min_segment_length. Finalmente, se asigna a cada punto una coordenada vertical acorde con la altura de referencia de vuelo y se construye un mensaje nav_msgs/Path con la trayectoria suavizada. La Figura 6.5 ilustra el efecto de este procedimiento sobre una ruta típica: la trayectoria poligonal generada por el RRT se representa junto a la curva resultante del suavizado por Chaikin. La B-spline conserva el corredor de espacio libre definido por la ruta original, pero reduce los ángulos de giro y distribuye la curvatura de forma más homogénea, lo que
CAPÍTULO 6. PERCEPCIÓN VISUAL PARA NAVEGACIÓN AUTÓNOMA
se traduce en referencias más adecuadas para el seguimiento por parte del cuadricóptero. Figura 6.5. Ruta poligonal generada por el planificador RRT y trayectoria suavizada obtenida a partir del algoritmo de Chaikin en el mismo entorno.
6.6.
Seguimiento de trayectorias y máquina de estados de misión
6.6.1.
Descripción general del seguidor de trayectorias El seguimiento de la trayectoria suavizada se realiza en un nodo dedicado de ROS, que actúa como interfaz entre la planificación en el espacio y el controlador del cuadricóptero. Este nodo recibe la trayectoria B-spline publicada en /dynamicNavigation/ bspline_trajectory y la odometría del vehículo, disponible en un tópico de tipo nav_msgs/Odometry; a partir de esta información genera consignas de posición y orientación en el marco map, que se publican en /CERLAB/quadcopter/setpoint_ pose.
Además de las referencias de pose, el seguidor gestiona los eventos de despegue y aterrizaje mediante los tópicos /CERLAB/quadcopter/takeoff y /CERLAB/quadcopter/land, activados cuando así lo indican los parámetros de configuración auto_takeoff y auto_land. Asimismo, el nodo mantiene un tópico de estado de
misión /dynamicNavigation/mission_state, en el que se publican etiquetas simbólicas que describen la fase actual del vuelo. Estas etiquetas son utilizadas por otros módulos, como el planificador local, para ajustar su comportamiento según si el sistema se encuentra en reposo, en ascenso, en seguimiento de trayectoria o en maniobra de aterrizaje.
6.6.2.
Máquina de estados de misión La lógica interna del seguidor se organiza como una máquina de estados finitos con cuatro fases principales de misión: idle, taking_off, following y landing. La Figura 6.6 resume esta estructura y las transiciones entre estados asociadas a los eventos de planificación y a las condiciones de vuelo.
En el estado idle no existe una trayectoria activa y el cuadricóptero permanece en reposo a la espera de recibir un mensaje nav_msgs/Path en el tópico de trayectorias suavizadas. Cuando se recibe una trayectoria válida y la opción de despegue automático auto_takeoff está habilitada, el seguidor emite un comando de despegue y transiciona a taking_off. Durante esta fase se monitoriza la altura estimada y, una vez alcanzado un umbral configurable, se asume que el vehículo ha completado el despegue y se pasa al estado following.
En following se recorre secuencialmente la lista de puntos de la B-spline, generando consignas de pose que combinan la posición del siguiente waypoint con una orientación de guiñada coherente con la dirección de avance. La distancia al punto de referencia se evalúa de forma continua y, cuando cae por debajo de una tolerancia predefinida, se avanza al siguiente waypoint. Al alcanzar el último punto de la trayectoria, y siempre que la opción auto_land esté activada, el seguidor inicia la secuencia de aterrizaje y cambia al estado landing, en el que se generan referencias de descenso vertical controlado hasta alcanzar una altura próxima al suelo. Completada esta maniobra, el sistema retorna al estado idle, quedando preparado para una nueva misión.
6.6.3.
Actualización en línea de trayectorias La integración con el planificador local requiere que el seguidor sea capaz de actualizar la trayectoria de referencia durante el vuelo sin interrumpir la misión ni forzar un nuevo despegue. Cuando el nodo se encuentra en el estado following y recibe una nueva trayectoria suavizada, se reemplaza la lista interna de puntos por la nueva B-spline y se selecciona un punto de inserción sobre esta curva, coherente con la posición actual
CAPÍTULO 6. PERCEPCIÓN VISUAL PARA NAVEGACIÓN AUTÓNOMA
F idle taking_off landing Inicio 'nueva trayectoria' 'altura alcanzada' 'comando de aterrizaje' 'misión completada' nueva trayectoria following 'recalcula trayectoria' Figura 6.6. Máquina de estados del seguidor de trayectorias. Se representan los estados idle, taking_off, following y landing, junto con las transiciones asociadas a la recepción de nuevas trayectorias y a las condiciones de vuelo. del cuadricóptero. Denotando por γnew(s) la trayectoria parametrizada por un parámetro de arco s y por pUAV la posición actual del vehículo, el índice de inserción se define como s∗= arg m´ın s≥0
γnew(s) −pUAV
,
(6.4)
es decir, el punto de la nueva trayectoria que minimiza la distancia al cuadricóptero. En la implementación discreta, esta búsqueda se realiza sobre el conjunto de puntos muestreados de la B-spline.
Para evitar movimientos hacia atrás indeseados, la selección puede restringirse a valores de s por encima de un umbral sm´ın asociado al progreso ya realizado a lo largo de la trayectoria anterior, de modo que solo se consideren puntos por delante de la posición actual. Una vez determinado el índice de inserción, el seguidor actualiza el waypoint activo y continúa el seguimiento desde ese punto de la nueva curva, manteniéndose en el estado following. Esta lógica de actualización en línea permite incorporar de forma natural las replanificaciones generadas por el planificador local sin reinicializar la misión ni introducir discontinuidades bruscas en el movimiento del cuadricóptero.
6.7.
Planificación local reactiva y costura de trayectorias
6.7.1.
Ventana local y detección de conflicto La planificación global sobre el mapa voxelizado proporciona una trayectoria que evita los obstáculos conocidos, pero no garantiza por sí sola una respuesta adecuada frente a cambios en el entorno o incertidumbres en el mapeo. La planificación se formula en el espacio cartesiano tridimensional definido por el mapa de ocupación voxelizado, construido a partir de la información perceptual del entorno. Para incorporar capacidad de reacción a escala local se introduce un planificador que opera sobre una ventana tridimensional centrada en la posición instantánea del cuadricóptero. Denotando por (x0, y0, z0) la posición actual del vehículo, la ventana de trabajo se define como W = [x0 −L, x0 + L] × [y0 −L, y0 + L] × [zm´ın, zm´ax],
(6.5)
donde L fija la extensión horizontal de la región considerada y [zm´ın, zm´ax] delimita el intervalo vertical relevante para la maniobra de evasión. A partir de esta definición se construye una instancia de VoxelGrid3D restringida a W, en la que se proyectan únicamente los puntos del mapa global que caen dentro de la ventana. Sobre esta grilla local se evalúa la seguridad de la trayectoria suavizada en un horizonte de distancia lookahead_distance. Para ello se identifica el punto de la B-spline más próximo a la posición actual del cuadricóptero y, a partir de dicho punto, se recorre la trayectoria acumulando longitud de arco hasta alcanzar ese máximo. En cada punto muestreado se estima la distancia mínima a las celdas ocupadas de la grilla local mediante una función dobs(p). Un punto de la trayectoria se considera peligroso cuando satisface dobs(p) < dsafe,
(6.6)
donde dsafe incorpora el tamaño efectivo del cuadricóptero, representado mediante una caja de colisión (parámetro collision_box), y un margen adicional de seguridad. En la implementación este criterio se evalúa mediante la función point_too_close_to_ obstacle, que explora un vecindario de vóxeles alrededor de la proyección del punto p.
El primer punto de la trayectoria que viola la condición (6.6) marca el inicio de un tramo conflictivo. A partir de ese índice se continúa avanzando sobre la misma curva hasta localizar el primer punto que vuelve a ser seguro, que se adopta como punto de reconexión Pconn sobre la ruta global. La situación se esquematiza en la Figura 6.7, donde la trayectoria original presenta un segmento que se aproxima en exceso a los obstáculos
CAPÍTULO 6. PERCEPCIÓN VISUAL PARA NAVEGACIÓN AUTÓNOMA
dentro de la ventana local y se señalan el intervalo conflictivo y el punto de reconexión seleccionado.
Figura 6.7. Esquema de planificación local en una ventana tridimensional alrededor del cuadricóptero.
6.7.2.
Replanificación local acotada y costura de trayectorias Delimitado el segmento conflictivo y determinado el punto de reconexión Pconn, el problema local se formula como una planificación acotada dentro de la ventana W. El punto de partida de la replanificación es la posición actual del cuadricóptero, proyectada a la celda libre más cercana en la grilla local si es necesario, y el objetivo es la posición de Pconn expresada en el mismo marco. Con estos datos se instancia un planificador RRTPlanner2p5D restringido a los límites de W, utilizando tiempos de timeout reducidos y parámetros de expansión adaptados al tamaño de la ventana. El resultado es un camino local que conecta la posición actual con Pconn evitando las celdas ocupadas que originaban el conflicto. Para integrarlo en la misión global se aplica un procedimiento de costura de trayectorias (trajectory stitching). Denotando por {γi} la secuencia discreta de puntos de la B-spline original, por istart el índice del primer punto peligroso y por jconn el índice de Pconn, la nueva trayectoria se construye concatenando: el prefijo seguro {γ0, . . . , γistart−1}, el camino local calculado por el RRT restringido, el sufijo {γjconn+1, . . . , γN} de la trayectoria original.
Este flujo se ilustra en la Figura 6.8, que muestra una secuencia típica de replanificación local y costura sobre el mapa voxelizado.
Figura 6.8. Secuencia de replanificación local y costura de trayectorias: (i) ruta global suavizada, (ii) detección del tramo conflictivo en la ventana local, (iii) generación de un camino local alternativo mediante RRT restringido y (iv) trayectoria resultante tras la costura.
La arquitectura descrita articula la percepción visual, la planificación jerárquica y el seguimiento de trayectorias en torno a un mapa de ocupación voxelizado. El cuadricóptero dispone de referencias suaves generadas por el planificador global y refinadas por el suavizador B-spline, mientras el seguidor de trayectorias y el planificador local reaccionan a la evolución del entorno mediante actualización en línea de rutas y maniobras de esquiva acotadas. Este entramado de componentes constituye el marco operativo sobre el que, en los análisis posteriores, se evaluará cuantitativamente el comportamiento del sistema en diferentes escenarios de navegación y su interacción con las estrategias de control desarrolladas previamente.
Capítulo 7 Análisis y resultados
7.1.
Resultados del modelo ARX MIMO (orden seis, retardo unitario) Con el propósito de contrastar el desarrollo teórico del Capítulo 3, se compara la salida medida del cuadricóptero con la salida predicha a un paso por el modelo ARX multivariable de orden seis con retardo unitario estimado fuera de línea mediante Mínimos Cuadrados. El ensayo de identificación se formula sobre el modelo linealizado alrededor del punto de operación de hover, entendido como la condición de equilibrio en torno a la cual se construye la aproximación lineal. Sin embargo, la prueba se desarrolla en lazo abierto y bajo señales de excitación aplicadas sin una ley de control reguladora que mantenga la dinámica acotada alrededor de ese equilibrio durante todo el intervalo de simulación. El regresor Φ se construye según la Subsección 3.1 y la ecuación (3.10), y los parámetros se obtienen resolviendo ˆθ = (Φ⊤Φ)−1Φ⊤Y (ecuaciones (3.13)–(3.16)), reorganizados luego en { ˆAℓ, ˆBℓ}6 ℓ=1 de acuerdo con (3.11). Cuando el acondicionamiento lo requiere, la solución se calcula de forma estable (QR/SVD) o con regularización εI, sin modificar la lógica del estimador.
Durante el ensayo de identificación, la planta se excita mediante señales binarias pseudoaleatorias en configuración MIMO, aplicadas de forma independiente en los cuatro canales de entrada, con el fin de garantizar persistencia de excitación. La formulación de estas señales, así como sus parámetros de muestreo, duración y amplitud, se presentó en la Subsección 3.3.1 y se ilustra en la Figura 3.1. A partir de esta excitación se registran las seis salidas (x, y, z, ϕ, θ, ψ) con un mismo período de muestreo. La evaluación se realiza sobre el mismo conjunto de datos utilizado para la estimación (in-sample) y, por claridad, las salidas se presentan en dos paneles:
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
Measured vs Predicted (y1–y3) Measured Predicted Measured Predicted Time [s] Measured Predicted [m] [m] [m] (a) y1, y2, y3 Measured vs Predicted (y4–y6) Measured Predicted Measured Predicted Time [s] Measured Predicted [rad] [rad] [rad] (b) y4, y5, y6 Figura 7.1. Superposición medida–predicha (predicción de un paso) para el modelo ARX MIMO de orden seis con retardo unitario, estimado por mínimos cuadrados fuera de línea. La Figura 7.1 muestra un seguimiento prácticamente punto a punto. En y1–y2 el modelo reproduce las curvaturas amplias del tramo final; en y3 mantiene el paralelismo de una rampa; y en y4–y6, de menor amplitud relativa, los desfasajes iniciales son leves y se atenúan con rapidez, lo que sugiere que el término autorregresivo compensa adecuadamente el efecto de entrada. Esta lectura es consistente con los indicadores: FIT ≈99,97 % (véase (7.3)), R2 ≃1 y RMSE/MAE casi constantes (≈3,14 y ≈2,52), lo que apunta a residuos de baja energía y homogéneos entre salidas.
Desde el punto de vista temporal, los pequeños errores se concentran en el arranque y en cambios de curvatura, en coherencia con una predicción un paso adelante: el uso de retardos recientes de y y u proporciona información inmediata que reduce el error en régimen y en trayectorias suaves. Como comprobaciones adicionales, se evalúan la blancura del residuo y su baja correlación cruzada con las entradas en rezagos relevantes; además, se verifica la estabilidad discreta de A(z) = I + P6 ℓ=1 ˆAℓz−ℓy la consistencia de la realización en espacio de estados (estabilizabilidad/detectabilidad).
Tabla 7.1. Métricas por canal i (Nv muestras).
Métrica Definición RMSEi q
1
Nv PNv k=1 ei[k]2 MAEi
1
Nv PNv k=1 |ei[k]| MSEi
1
Nv PNv k=1 ei[k]2 = RMSE 2 i R2 i
1 −
PNv k=1 ei[k]2 PNv k=1(yi[k] −¯yi)2
7.1.1.
Métricas de evaluación y validación El conjunto considerado tiene Ts = 0,01 s y N = 6001 muestras; tras el alineamiento por retardos, Nv = 5994. Las métricas por canal se definen a partir del residuo ei[k] = yi[k]−ˆyi[k] y del promedio ¯yi =
1
Nv PNv k=1 yi[k]. La Tabla 7.1 resume las expresiones, y el FIT utilizado como referencia se define en (7.3). Las métricas globales se calculan como promedios simples entre canales (7.4).
R2 i = 1 − P k(yi[k] −ˆyi[k])2 P k(yi[k] −¯yi)2 ,
(7.1)
MSEi = 1 N Nv X k=1 (yi[k]−ˆyi[k])2, RMSEi = q MSEi, MAEi = 1 N Nv X k=1 |yi[k]−ˆyi[k]|,
(7.2)
FITi( %) = 100
1 −∥yi −ˆyi∥2
∥yi −¯yi∥2 !
,
(7.3)
RMSEglob = 1 p p X i=1 RMSEi, MAEglob = 1 p p X i=1 MAEi, R2 glob = 1 p p X i=1 R2 i , FITglob = 1 p p X i=1 FITi.
(7.4)
Los indicadores confirman un ajuste in-sample muy elevado (FIT ≈99,9 %, R2 ≃
1) y residuos de baja amplitud, coherentes con la superposición medida–predicha de la
Figura 7.1. Las trazas de error no muestran patrones evidentes y su inspección visual sugiere blancura; adicionalmente, A(z) cumple estabilidad en el rango de interés y la realización en espacio de estados es consistente con los criterios de estabilizabilidad y detectabilidad.
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
Tabla 7.2. Métricas de ajuste (predicción de un paso) para el modelo ARX MIMO de orden seis, retardo unitario.
Canal
RMSE
MAE
R2
FIT ( %)
MSE
Global
3.143
2.517
1.0000
99.97
9,881 × 100 y1
3.143
2.517
1.0000
100.00
9,881 × 100 y2
3.143
2.517
1.0000
100.00
9,881 × 100 y3
3.143
2.517
1.0000
99.94
9,881 × 100 y4
3.143
2.517
1.0000
99.98
9,881 × 100 y5
3.143
2.517
1.0000
99.94
9,881 × 100 y6
3.143
2.517
1.0000
99.96
9,881 × 100 Cifras redondeadas; Nv = 5994 muestras válidas.
0
10
20
30
40
50
60
-6
-4
-2
0
y1
106
Measured vs Estimated Outputs (y1–y3) Measured Estimated
0
10
20
30
40
50
60
-4
-2
0
y2
106
Measured Estimated
0
10
20
30
40
50
60
Time [s]
0
50
100
y3 Measured Estimated [m] [m] [m] (a) y1, y2, y3
0
10
20
30
40
50
60
0
200
400
600
y4 Measured vs Estimated Outputs (y4–y6) Measured Estimated
0
10
20
30
40
50
60
-1000
-500
0
y5 Measured Estimated
0
10
20
30
40
50
60
Time [s]
-400
-200
0
y6 Measured Estimated [m] [m] [m] (b) y4, y5, y6 Figura 7.2. Estimación en línea (RLS) para ARX MIMO de orden seis y retardo unitario — predicción a un paso.
7.2.
Resultados de la estimación en línea mediante RLS Con el fin de evaluar el desempeño en línea del estimador presentado en la Sección 3.3.1, se contrasta la salida medida del cuadricóptero con la salida estimada a un paso por RLS aplicado al modelo ARX multivariable de orden seis con retardo unitario. La actualización paramétrica sigue (3.19)–(3.20), mientras que el regresor se forma en tiempo real según (3.9)–(3.10). La simulación se ejecuta con Ts = 0,01 s durante N = 6001 muestras, bajo entradas PRBS independientes por canal.
Las subfiguras de la Figura 7.2 agrupan las seis salidas para favorecer la lectura. En y1–y2 el estimador reproduce con precisión las variaciones de mayor escala; en y3 se conserva el paralelismo propio de trayectorias tipo rampa; y en y4–y6, de menor amplitud
-2
0
2
e1(t) Absolute Estimation Error (e1–e3)
0
10
20
30
40
50
60
-2
0
2
4
e2(t)
0
10
20
30
40
50
60
Time [s]
-2
0
2
e3(t) (a) e1, e2, e3.
0
10
20
30
40
50
60
-2
0
2
e4(t) Absolute Estimation Error (e4–e6)
0
10
20
30
40
50
60
-2
0
2
e5(t)
0
10
20
30
40
50
60
Time [s]
-2
0
2
e6(t) (b) e4, e5, e6.
Figura 7.3. Errores absolutos de estimación por salida en el esquema RLS–ARX MIMO.
0
10
20
30
40
50
60
-10
0
10
erel,1 [%] Relative Estimation Error (erel,1–erel,3)
0
10
20
30
40
50
60
-10
0
10
erel,2 [%]
0
10
20
30
40
50
60
Time [s]
-10
0
10
erel,3 [%] (a) erel,1, erel,2, erel,3 ( %).
0
10
20
30
40
50
60
-10
0
10
erel,4 [%] Relative Estimation Error (erel,4–erel,6)
0
10
20
30
40
50
60
-10
0
10
erel,5 [%]
0
10
20
30
40
50
60
Time [s]
-10
0
10
erel,6 [%] (b) erel,4, erel,5, erel,6 ( %).
Figura 7.4. Errores relativos de estimación por salida ( %) en el esquema RLS–ARX MI- MO.
relativa, los desfasajes iniciales son reducidos y se atenúan con rapidez, en consonancia con la contribución autorregresiva del predictor.
En las subfiguras de la Figura 7.3 se observan los errores absolutos ei(t) separados por grupos de salidas. Tras el arranque, los ei(t) se concentran cerca de cero, sin patrones persistentes, y se comportan como ruido de varianza acotada. En e1–e3 aparecen fluctuaciones algo mayores alrededor de t ≈35–45 s, que coinciden con cambios de curvatura en las salidas, mientras que e4–e6 permanecen más compactos y estables en todo el horizonte. En conjunto, la energía del residuo es baja y consistente con el ajuste cuantitativo reportado por RMSE y MAE.
La Figura 7.4 muestra los errores relativos erel,i(t), normalizados por la energía de cada salida. Esta normalización colapsa las diferencias de escala y revela bandas estrechas en torno a 0 % para todos los canales (típicamente dentro de ± 1 % a 2 %). Tras el período
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
Tabla 7.3. Métricas de ajuste (predicción a un paso) — RLS en línea para ARX MIMO de orden seis.
Canal
RMSE
MAE
R2
FIT ( %)
MSE
Global
0.38539
0.30407
1.0000
99.76
1,50292 × 10−1 y1
0.43783
0.34391
1.0000
100.00
1,91699 × 10−1 y2
0.45135
0.34949
1.0000
100.00
2,03715 × 10−1 y3
0.35567
0.28270
0.9999
99.22
1,26500 × 10−1 y4
0.35581
0.28264
1.0000
99.86
1,26602 × 10−1 y5
0.35597
0.28279
1.0000
99.92
1,26713 × 10−1 y6
0.35570
0.28289
1.0000
99.58
1,26521 × 10−1 Cifras reportadas para N=6001 y predicción a un paso.
de calentamiento, no se aprecian derivas ni oscilaciones sistemáticas, lo que sugiere una sintonía adecuada del factor de olvido y correcciones paramétricas oportunas. El error absoluto ei(t) cuantifica la desviación del modelo en unidades físicas de la salida, mientras que el error relativo erel,i(t) expresa la magnitud de esa desviación en proporción a la propia señal. En este caso, los errores absolutos se mantienen bajos y las bandas relativas son estrechas; ambas lecturas concuerdan con métricas altas de desempeño (FIT y R2) y respaldan la fidelidad del predictor a un paso en el régimen analizado. Esta lectura visual se complementa con las métricas definidas en (7.1), sintetizadas en la Tabla 7.3.
Los valores de la Tabla 7.3 evidencian un ajuste in-sample muy elevado. En particular, FITglob = 99,76 % y R2 ≃1 indican que la energía del residuo es sustancialmente menor que la variabilidad de las señales ((7.1)–(7.3)). La proximidad entre RMSE y MAE sugiere errores acotados y relativamente homogéneos; además, la relación MSE ≈RMSE2 confirma la consistencia numérica ((7.2)). Por canales, y1 y y2 alcanzan un FIT prácticamente del 100 %; y3 presenta el valor más bajo (99,22 %), compatible con ligeras discrepancias en tramos de mayor curvatura; mientras que y4–y6 mantienen desempeños cercanos o superiores al 99,5 % a lo largo del horizonte.
7.3.
Resultados del controlador nominal LQI Se evalúa el seguimiento en [x, y, z, ψ] con acción integral (LQI) en tiempo discreto, bajo Ts = 0,01 s y referencias escalonadas en instantes distintos. Las Figuras 7.6–7.8 muestran respuesta, entradas y errores; la Tabla 7.5 resume métricas canal a canal.
r[k] = xd yd zd ψd + − e[k]
1
z − 1 Integrador discreto Ki + − u[k] B + +
1
z − 1 x[k] C y[k] = x y z ψ A K y[k] Realimentación de estados r[k]: referencia e[k] : error de seguimiento Ki : ganancia integral K : ganancia óptima LQR x[k], y[k], u[k]: estado, salida y control Planta discreta Figura 7.5. Esquema resumido del lazo LQI nominal utilizado en la evaluación de resultados.
Tabla 7.4. Variables empleadas en el lazo LQI nominal.
Tipo Variables Descripción Referencia r[k] = [xd, yd, zd, ψd]⊤ Consignas de posición y guiñada Salida y[k] = [x, y, z, ψ]⊤ Salidas reguladas y medidas Error e[k] = r[k] −y[k] Error de seguimiento Control u[k] = [T, τx, τy, τz]⊤ Empuje total y torques aplicados Ganancias Kx, Ki Ganancias nominales del LQI
7.3.1.
Estructura del lazo y variables empleadas La Figura 7.5 resume las variables empleadas en la evaluación del controlador nominal LQI, cuyo diseño se desarrolló en la Subsección 4.2.2. El lazo regula las salidas y[k] = [x, y, z, ψ]⊤a partir de la referencia r[k] = [xd, yd, zd, ψd]⊤, el error de seguimiento e[k] y el estado integral ξ[k], generando la señal de control u[k] = [T, τx, τy, τz]⊤. En la Figura 7.6 las salidas convergen sin error estacionario, conforme a la acción integral. El canal x alcanza rápidamente su referencia con leve sobreimpulso; y presenta un sobreimpulso elevado alrededor del escalón y se amortigua sin oscilación sostenida; z es más lento por la jerarquía actitud–posición y el instante tardío del escalón; ψ converge de forma suave, aunque con el mayor tiempo de establecimiento. Las entradas (Figura 7.7) muestran picos bien localizados al inicio de cada transición y un retorno rápido al régimen, manteniéndose dentro de márgenes razonables. Los errores (Figura 7.8) decaen monótonamente tras cada referencia, con colas cortas y sin derivas.
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
0
2
4
6
8
10
12
14
16
18
0
1
2
x Tracking x, y, z, Output Reference
0
2
4
6
8
10
12
14
16
18
-2
-1
0
y
0
2
4
6
8
10
12
14
16
18
0
1
z
0
2
4
6
8
10
12
14
16
18
Time [s]
0
0.5
1
psi [m] [m] [m] [rad] Figura 7.6. Seguimiento de referencias en x, y, z, ψ con LQI.
7.3.2.
Métricas de evaluación y validación La Tabla 7.5 resume el desempeño del LQI nominal por canal. En términos de error, ψ presenta los valores más bajos de IAE, ISE y RMSE, mientras que x y y muestran errores globales similares. El canal z registra el mayor ITAE debido al instante tardío de su cambio de referencia, lo que incrementa la penalización temporal del error. El sobreimpulso porcentual OS (Overshoot) es reducido en x, z y ψ, pero elevado en y, asociado al acoplamiento entre la dinámica lateral y la actitud del vehículo. Finalmente, los tiempos de establecimiento test indican una respuesta más rápida en x y una convergencia más lenta en ψ. En conjunto, las métricas confirman un seguimiento estable y sin error estacionario, coherente con la acción integral del LQI.
-1
0
1
2
T Control Inputs
0
2
4
6
8
10
12
14
16
18
-0.2
0
0.2
x
0
2
4
6
8
10
12
14
16
18
-0.2
0
0.2
0.4
y
0
2
4
6
8
10
12
14
16
18
Time [s]
0
0.1
0.2
z [N] [N m] [N m] [N m] Figura 7.7. Entradas de control u = [T, τx, τy, τz]⊤durante el seguimiento. Tabla 7.5. Métricas por canal para el LQI nominal.
Canal
IAE
ISE
ITAE
RMSE
OS ( %) test (s) x
1.890
2.945
4.807
0.3837
2.9
3.96
y
1.870
2.906
8.488
0.3812
100.0
5.95
z
1.725
1.897
11.603
0.3080
4.0
8.37
ψ
0.348
0.195
2.884
0.0987
2.4
9.80
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
0
2
4
6
8
10
12
14
16
18
0
1
2
ex(t) Tracking Errors
0
2
4
6
8
10
12
14
16
18
-2
-1
0
ey(t)
0
2
4
6
8
10
12
14
16
18
0
0.5
1
ez(t)
0
2
4
6
8
10
12
14
16
18
Time [s]
0
0.2
0.4
0.6
eA(t) Figura 7.8. Errores de seguimiento ex, ey, ez, eψ.
7.4.
Resultados del Control Adaptativo Óptimo — STR Indirecto Esta sección recoge los principales resultados del esquema STR indirecto basado en LQI nominal seguro, identificación RLS–ARX MIMO y rediseño adaptativo del LQI con ganancia de referencia. Se emplea una referencia multi–step (escalonada) sobre las salidas reguladas [x, y, z, ψ], con horizonte de simulación de 40 s y cambio paramétrico en t = 20 s. Durante el arranque, el LQI nominal físico de 12 estados actúa como modo seguro, mientras que el identificador en línea habilita la transición inicial hacia el modo adaptativo, los reajustes supervisados de ganancias y, cuando es necesario, el retorno seguro al modo nominal.
7.4.1.
Estructura del esquema STR y variables empleadas La Figura 7.9 resume la estructura general del esquema STR indirecto utilizada para la evaluación. En ella se distinguen las referencias de seguimiento, las salidas reguladas, las entradas de control y las variables asociadas al proceso de identificación y reajuste supervisado.
7.4.2.
Métricas de evaluación y validación A continuación se sintetiza el desempeño del STR indirecto. Primero se reportan las métricas de seguimiento por canal y, después, un conjunto de indicadores compactos que caracterizan el lazo LQI y la calidad del modelo RLS.
La Tabla 7.7 resume el desempeño del STR indirecto bajo la referencia multi–step. Las métricas por canal permiten identificar la contribución de cada salida al error total, mientras que la fila global resume el comportamiento conjunto del sistema. Se observa que los mayores aportes al error provienen de los canales x y y, asociados a los cambios de referencia en posición horizontal, mientras que z y ψ mantienen errores más contenidos. La energía de control Eu permanece en un nivel reducido, lo que evidencia que la adaptación no requiere acciones excesivamente agresivas.
La Tabla 7.8 resume los eventos principales del supervisor. La transición inicial hacia el modo adaptativo se realiza en t = 27,49 s, una vez que el identificador cumple los criterios de calidad definidos. A partir de ese instante, el controlador puede reajustar sus ganancias de forma supervisada: se aceptan 58 reajustes y se rechazan 15 intentos, lo que confirma que la adaptación no se aplica de manera indiscriminada, sino condicionada por
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
LAZO DE ADAPTACIÓN (STR INDIRECTO)
LAZO DE CONTROL PRINCIPAL (LQI)
Diseño del Controlador (cálculo de ganancias) Resolución DARE ⇒ Kx[k], Ki[k] Estimación RLS (Identificación en línea) Modelo ARX MIMO y[k] = Σ Ai y[k-i] + Σ Bj u[k-j] + e[k] θ̂[k] Parámetros estimados (Âi, B̂j) r[k] = ⎡ xd ⎤ ⎢ yd ⎥ ⎢ zd ⎥ ⎣ ψd ⎦ + − e[k]
1
z − 1 Integrador discreto ξ[k] Estado integral Controlador LQI (con ganancias K[k]) u[k] = −Kx[k] x[k] − Ki[k] ξ[k] u[k] ⎡ T ⎤ ⎢ τx ⎥ ⎢ τy ⎥ ⎣ τz ⎦ Planta Quadcopter y[k] ⎡ x ⎤ ⎢ y ⎥ ⎢ z ⎥ ⎣ ψ ⎦ u[k] y[k] Kx[k] Ki[k] Ganancias actualizadas por STR Lazo de control principal (LQI) Lazo de adaptación (STR indirecto) r[k]: referencia e[k]: error de seguimiento ξ[k]: estado integral u[k]: señal de control y[k]: salidas medidas θ̂[k]: parámetros estimados Kx[k], Ki[k]: ganancias LQI K[k]: ganancias actualizadas − Figura 7.9. Esquema resumido del STR indirecto utilizado en la evaluación de resultados. Tabla 7.6. Variables empleadas en el esquema LQI–STR indirecto. Tipo Variables Descripción Referencia r[k] = [xd, yd, zd, ψd]⊤ Consignas de posición y guiñada Salida medida y[k] = [x, y, z, ψ]⊤ Salidas reguladas del cuadricóptero Error e[k] = r[k] −y[k] Error de seguimiento Control u[k] = [T, τx, τy, τz]⊤ Empuje total y torques aplicados Identificación ˆθ[k], ˆAi, ˆBj Parámetros estimados del modelo ARX MIMO Ganancias Kx[k], Ki[k], Kr[k] Realimentación, acción integral y ganancia de referencia validaciones del modelo identificado, restricciones de actuación y ventanas de seguridad. El retorno seguro al modo nominal registrado evidencia que el LQI físico de 12 estados permanece disponible como respaldo cuando se detectan condiciones de riesgo.
7.4.3.
Referencias aplicadas La Figura 7.10 muestra los cambios programados de referencia que estructuran la evaluación. Estos escalones permiten evaluar la respuesta del controlador ante consignas sucesivas, sin saturar sostenidamente los actuadores ni el integrador.
Tabla 7.7. Métricas de desempeño globales y por canal con STR indirecto. Salida
IAE
ISE
ITAE
RMSE
Eu x
3.6286
3.2385
49.686
0.2845
– y
4.5330
7.9309
89.897
0.4453
– z
2.0634
1.0357
18.322
0.1609
– ψ
1.5872
1.6725
30.681
0.2045
– Global / Total
11.812
13.878
188.59
0.2738
2.3068
Tabla 7.8. Indicadores de supervisión del STR indirecto.
Indicador Valor Transición inicial al modo adaptativo Sí tadapt [s]
27.49
Reajustes aceptados
58
Reajustes rechazados
15
Retornos seguros al modo nominal
1
Saturación promedio [ %]
0.025
Máx. |ϕ| [deg]
42.533
Máx. |θ| [deg]
20.564
7.4.4.
Seguimiento de posición y guiñada La Figura 7.11 muestra el seguimiento compacto de posición y guiñada. Las salidas siguen la referencia multi–step con transitorios acotados y error final prácticamente nulo. La mayor exigencia se observa en los canales x y y, donde los cambios de referencia demandan mayor acción de actitud, mientras que z y ψ presentan una respuesta más contenida.
7.4.5.
Arranque seguro con LQI 12×12 La Figura 7.12 evidencia el comportamiento durante la fase de arranque con el LQI nominal físico de 12 estados. Se muestra hasta t = 27,5 s, ya que la transición inicial hacia el modo adaptativo ocurre en tadapt = 27,49 s, según la Tabla 7.8. Este modo seguro permanece disponible durante toda la simulación como respaldo ante modelos no aceptados, saturación reciente o condiciones de riesgo.
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
0
5
10
15
20
25
0
0.5
1
1.5
2
xr [m] Referencia xr
0
5
10
15
20
25
-2
-1.5
-1
-0.5
0
0.5
1
yr [m] Referencia yr
0
5
10
15
20
25
t [s]
1
1.1
1.2
1.3
1.4
zr [m] Referencia zr
0
5
10
15
20
25
t [s]
-60
-40
-20
0
20
40
60
Ar [°] Referencia Ar Figura 7.10. Referencias seleccionadas para x, y, z y ψ (multi–step).
7.4.6.
Evolución completa de las seis salidas En la Figura 7.13 se aprecia la coordinación actitud–posición: ϕ y θ reaccionan como lazos internos rápidos que facilitan el seguimiento de x y y; las oscilaciones se atenúan sin efectos persistentes, y ψ acompasa su referencia con sobrepasos moderados.
7.4.7.
Errores de seguimiento La Figura 7.14 muestra que, tras cada escalón, los errores decaen rápidamente hacia cero y se mantienen dentro de bandas estrechas. El mayor aporte al IAE/ISE proviene de las transiciones en y; z y ψ exhiben errores de menor energía y prácticamente nulo sesgo final.
7.4.8.
Lectura conjunta de métricas y supervisión Las figuras y tablas muestran que el STR indirecto mantiene seguimiento acotado en los canales regulados y error final prácticamente nulo. La acción integral del LQI permite
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
Posición [m] Seguimiento de posición (X,Y,Z) x ref_x y ref_y z ref_z
0
5
10
15
20
25
30
t [s]
-40
-20
0
20
40
60
Yaw [°] A ref_A Figura 7.11. Seguimiento de x, y, z (arriba) y ψ (abajo). Se sobreponen referencias (líneas punteadas).
eliminar el sesgo estacionario, mientras que la ganancia de referencia contribuye a mejorar la respuesta frente a cambios de consigna. A su vez, el supervisor evita aplicar reajustes cuando el modelo identificado o las condiciones de actuación no son favorables. En conjunto, el esquema combina tres elementos: un LQI nominal físico como modo seguro, un identificador RLS–ARX MIMO que actualiza el modelo de los canales regulados y un rediseño adaptativo del LQI con ganancia de referencia. Esta combinación permite obtener un desempeño competitivo frente a controladores de ganancias fijas, conservando mecanismos de protección como bloqueo por saturación reciente, validación de reajustes supervisados y retorno seguro al modo nominal.
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
0
5
10
15
20
25
30
-0.5
0
0.5
1
1.5
2
2.5
x [m] x vs xr (0--27.5 s) x xr
0
5
10
15
20
25
30
-2
-1.5
-1
-0.5
0
0.5
1
y [m] y vs yr (0--27.5 s) y yr
0
5
10
15
20
25
30
t [s]
0
0.5
1
1.5
z [m] z vs zr (0--27.5 s) z zr
0
5
10
15
20
25
30
t [s]
-60
-40
-20
0
20
40
60
80
A [°] A vs Ar (0--27.5 s) A Ar Figura 7.12. Arranque seguro con LQI 12×12: comparación salida–referencia por canal. Tabla 7.9. Comparación global de desempeño entre controladores. Controlador
IAE
ISE
ITAE
RMSE
Eu
STR–LQI
11.812
13.878
188.59
0.2738
2.3068
LQR–AG+Kr
11.591
13.994
195.42
0.2804
28.498
LQG+Kr
16.168
16.150
242.00
0.3261
0.1900
7.4.9.
Comparación global con controladores de referencia Con el fin de contextualizar el desempeño del esquema propuesto, la Tabla 7.9 resume una comparación global frente a dos controladores de referencia desarrollados en los apéndices: el LQR sintonizado mediante algoritmo genético con ganancia de referencia Kr (Apéndice G, Tabla G.3) y el LQG con Kr (Apéndice H, Tabla H.2). La comparación se realiza a partir de métricas integrales de error y energía de control acumulada Eu, considerando que valores menores indican mejor desempeño para cada indicador. La Tabla 7.9 muestra que el STR–LQI obtiene el mejor desempeño en ISE, ITAE y RMSE, lo que indica una reducción de la energía del error, menor penalización temporal y mejor ajuste global de seguimiento. El LQR–AG+Kr presenta el menor IAE, aunque
Evolución de todas las salidas (STR LQI--RLS--ARX)
0
5
10
15
20
25
30
35
40
0
0.5
1
1.5
2
x [m] x [m] x ref_x
0
5
10
15
20
25
30
35
40
-2
-1
0
1
y [m] y [m] y ref_y
0
5
10
15
20
25
30
35
40
0
0.5
1
1.5
z [m] z [m] z ref_z
0
5
10
15
20
25
30
35
40
-40
-20
0
20
? [°] ? [°] ?
0
5
10
15
20
25
30
35
40
t [s]
-20
-10
0
10
20
3 [°]
3 [°]
3
0
5
10
15
20
25
30
35
40
t [s]
-50
0
50
A [°] A [°] A ref_A Figura 7.13. Evolución de las 6 salidas del cuadricóptero (posición y actitud). con una energía de control considerablemente mayor. Por su parte, el LQG+Kr alcanza el menor Eu, pero a costa de mayores errores acumulados. En conjunto, el STR–LQI ofrece el desempeño global más equilibrado, al combinar precisión de seguimiento, esfuerzo de control moderado y capacidad de reajuste supervisado ante variaciones del modelo.
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
0
5
10
15
20
25
30
35
40
-2
-1
0
1
ex [m] Errores de seguimiento (x,y,z,A) ex
0
5
10
15
20
25
30
35
40
-1
0
1
2
3
ey [m] ey
0
5
10
15
20
25
30
35
40
-0.5
0
0.5
1
ez [m] ez
0
5
10
15
20
25
30
35
40
t [s]
-100
-50
0
50
100
eA [°] eA Figura 7.14. Errores ex, ey, ez y eψ a lo largo de la prueba.
7.5.
Resultados de navegación y planificación autónoma En esta sección se evalúa cuantitativamente el comportamiento de la arquitectura de navegación basada en mapas de ocupación voxelizados descrita en el Capítulo 6, dentro del entorno de simulación ROS–Gazebo/RViz presentado en la Sección. 6.2. El análisis se organiza en tres ejes complementarios: en primer lugar, la calidad geométrica de las trayectorias generadas por el planificador global y el suavizador B-spline; en segundo lugar, el desempeño de seguimiento y control del cuadricóptero frente a dichas trayectorias; y, finalmente, el rendimiento computacional de la planificación global y local a partir de las estadísticas recogidas durante las ejecuciones.
La Figura 7.15 presenta una secuencia cronológica representativa de la navegación autónoma en el entorno dinámico simulado. En cada panel se muestra, a la izquierda, la visualización en RViz del mapa de ocupación voxelizado y la trayectoria activa del cuadricóptero; y, a la derecha, la ejecución simultánea en Gazebo con obstáculos estáticos y
peatones móviles. La secuencia permite observar la detección de conflictos sobre la trayectoria activa, la solicitud de replaneo local, la actualización de la trayectoria seguida y la continuación de la misión una vez superado el evento. Esta evidencia visual complementa la interpretación cualitativa del comportamiento de evasión de obstáculos del sistema y respalda las métricas cuantitativas resumidas en la Tabla 7.10. i ii iii iv Figura 7.15. Secuencia cronológica de navegación autónoma en un entorno interior dinámico. (i) Detección de conflicto sobre la trayectoria activa y solicitud de replaneo local. (ii) Recepción de una trayectoria alternativa para evasión del obstáculo. (iii) Seguimiento de la trayectoria actualizada tras el replaneo. (iv) Continuación de la misión hacia el objetivo de navegación.
Para cuantificar el desempeño de la arquitectura de navegación en el Escenario 1, se procesaron n = 3 recorridos independientes (simulación) registrados en rosbag, a partir de los cuales se extrajeron trayectorias planificadas (RRT y suavizado B-spline), odometría y referencias. Con estos datos se calcularon métricas geométricas (longitud, suavidad y curvatura), de seguimiento (RMSE de posición y guiñada) y de cómputo (tiempos de planificación global/local y eventos de replanificación). En este marco, los estadísticos reportados se interpretan como una caracterización descriptiva del comportamiento observado en el escenario evaluado. La Tabla 7.10 resume cada métrica mediante (i) media ± desviación estándar (σ) y (ii) mediana con [Q1, Q3] y rango (m´ın – m´ax), proporcionando una síntesis compacta de la tendencia central y la dispersión entre recorridos en ROS–Gazebo/RViz.
CAPÍTULO 7. ANÁLISIS Y RESULTADOS
Para evitar ambigüedades en la interpretación, los factores de mejora JRRT/JBS y κRRT/κBS se computan de forma individual en cada recorrido como cocientes entre la trayectoria RRT y su correspondiente trayectoria suavizada, y posteriormente se agregan con los estadísticos reportados en la Tabla 7.10. De manera análoga, los tiempos ¯Tglobal y ¯Tlocal corresponden a promedios por llamada (tiempos medidos en cada invocación del planificador registrada en el rosbag), y no a tiempos acumulados a lo largo de la misión. Tabla 7.10. Métricas de navegación (Escenario 1, n = 3 recorridos). Bloque Métrica Símbolo Media±σ Mediana [Q1,Q3] min–max Calidad geométrica Longitud de trayectoria suavizada [m]
LBS
4,762 ± 3,490
5.118
[3.113, 6.589] 1.107–
8.059
Índice de suavidad (B-spline) [−]
JBS
0,039 ± 0,021
0.039
[0.029, 0.050] 0.018–
0.060
Factor de mejora en suavidad [−]
JRRT/JBS 7,665 ± 5,485
6.093
[4.615, 9.928] 3.137–
13.764
Curvatura máxima (B-spline) [rad] κBS 0,515 ± 0,177
0.459
[0.416, 0.586] 0.373–
0.713
Factor de reducción de curvatura [−] κRRT/κBS 2,009 ± 0,127
2.061
[1.962, 2.082] 1.864–
2.102
Seguimiento y control RMSE de posición [m] RMSEpos 0,634 ± 0,094
0.623
[0.585, 0.678] 0.547–
0.733
Error relativo de posición [ % de
LBS]
epos,rel 28,396 ±
32,878
12.171
[9.478, 39.202] 6.785–
66.233
RMSE de guiñada [rad] RMSE(rad) ψ 0,472 ± 0,123
0.421
[0.402, 0.516] 0.383–
0.612
RMSE de guiñada [◦]
RMSE(◦)
ψ 27,031±7,035
24.109
[23.018, 29.582] 21.928–
35.056
Rendimiento de planificación (global y local) Tiempo medio de planificación global [s] ¯Tglobal 15,037±3,695
15.129
[13.213, 16.908] 11.297–
18.686
Tiempo medio de planificación local [s] ¯Tlocal 0,388 ± 0,054
0.418
[0.372, 0.419] 0.326–
0.419
Relación global/local [−] ¯Tglobal/ ¯Tlocal 39,358 ±
10,781
44.680
[35.815, 45.561] 26.951–
46.442
Tiempo mínimo de replanificación local [ms] T m´ın local 1,323 ± 0,061
1.337
[1.297, 1.356] 1.257–
1.376
En calidad geométrica, la Tabla 7.10 indica que el suavizado B-spline mejora de forma consistente la regularidad de la trayectoria frente a la ruta base del RRT: JRRT/JBS
presenta una mediana de 6,093 y κRRT/κBS se mantiene cercano a 2. La variabilidad observada en JRRT/JBS (rango 3,137–13,764) sugiere que la ganancia en suavidad depende de la geometría particular de cada misión, mientras que la reducción de curvatura muestra una dispersión acotada, compatible con un efecto de suavizado estable. En seguimiento y control, según la Tabla 7.10, la precisión de posición se resume con RMSEpos, cuya mediana es 0,623 m (rango 0,547–0,733 m). El error relativo epos,rel presenta mayor dispersión porque normaliza por la longitud de misión LBS, la cual varía entre recorridos; por ello se interpreta como métrica complementaria al RMSE absoluto. Para guiñada, RMSEψ exhibe una mediana de 0,421 rad (24,1◦), coherente con una validación preliminar en simulación donde la orientación se regula para mantener alineación hacia waypoints consecutivos.
En rendimiento computacional, la Tabla 7.10 refleja la diferencia de escalas entre capas: la planificación global presenta mayores tiempos medios (mediana ¯Tglobal = 15,129 s; rango 11,297–18,686 s), mientras que la planificación local permanece en el orden de décimas de segundo (mediana ¯Tlocal = 0,418 s; rango 0,326–0,419 s). El tiempo mínimo de replanificación local se mantiene en el orden de 1,3 ms, lo que sugiere una capacidad de reacción rápida ante cambios locales del mapa voxelizado en el entorno de simulación.
Capítulo 8 Conclusiones y trabajos futuros La investigación abordó la navegación autónoma de cuadricópteros en entornos tridimensionales parcialmente conocidos, con incertidumbre paramétrica y restricciones físicas y computacionales, mediante un sistema integrado que enlaza modelado dinámico, control óptimo en tiempo discreto con adaptación en línea y percepción basada en mapas de ocupación voxelizados. Esta integración se implementó sobre una plataforma embebida con simulación en ROS–Gazebo, manteniendo trazabilidad entre las hipótesis de modelado, las decisiones de diseño y el comportamiento observado en los escenarios de navegación. En conjunto, el trabajo consolida un flujo de navegación autónoma extremo a extremo, formulado para implementación embebida y documentado de manera reproducible.
En el ámbito del control, el avance principal consiste en la formulación y validación de un regulador LQI discreto MIMO sujeto a restricciones de muestreo y saturación, junto con su extensión adaptativa mediante STR indirecto. Los resultados en lazo cerrado muestran errores de seguimiento acotados y esfuerzos de control compatibles con los actuadores, mientras que la capa adaptativa mitiga variaciones paramétricas sin comprometer la estabilidad interna. Además, la comparación frente a LQR–AG+Kr y LQG+Kr evidencia que el STR–LQI ofrece el mejor desempeño global, al equilibrar precisión de seguimiento, esfuerzo de control moderado y adaptación supervisada frente a variaciones del modelo. De este modo, la investigación contribuye al control de vehículos aéreos no tripulados mediante un esquema LQI+STR discreto, supervisado y orientado a arquitecturas de cómputo embarcado.
En la dimensión de percepción y planificación, la arquitectura implementada consolida una navegación jerárquica apoyada en un mapa de ocupación voxelizado compartido por los módulos de decisión, que integra planificación global, suavizado geométrico de trayectorias y un planificador local reactivo, coordinados mediante una máquina de estados
CAPÍTULO 8. CONCLUSIONES Y TRABAJOS FUTUROS
de misión. Esta arquitectura permite generar trayectorias libres de colisión y compatibles con la dinámica del vehículo, así como ajustar localmente la ruta ante violaciones del margen de seguridad, manteniendo la continuidad de la misión. En consecuencia, el sistema aporta una solución de percepción y planificación diseñada para respetar las restricciones dinámicas del cuadricóptero y facilitar su integración con controladores de tipo LQI. La evaluación conjunta del sistema muestra que la interacción entre percepción voxelizada, planificación jerárquica y control LQI/STR es consistente con las restricciones físicas y computacionales adoptadas. El cuadricóptero sigue trayectorias suavizadas respetando márgenes de seguridad frente a obstáculos y esfuerzos de control compatibles con los actuadores, mientras que la capa adaptativa introduce robustez frente a variaciones de modelo sin interferir con el ciclo nominal ni con las salvaguardas de saturación. Frente a enfoques que analizan por separado el control de actitud y posición o que formulan la navegación sobre modelos cinemáticos simplificados, el trabajo sitúa un puente operativo entre ambas perspectivas mediante un banco de pruebas reproducible que hace explícitas las decisiones de diseño y su impacto en el comportamiento del sistema. En síntesis, el sistema desarrollado demuestra que la combinación de control óptimo con adaptación indirecta y navegación basada en mapas voxelizados permite abordar de forma integrada la regulación, el seguimiento de trayectorias y la evasión de obstáculos en entornos tridimensionales controlados, bajo restricciones de cómputo y de actuación comparables a las de plataformas reales. Esta demostración constituye la contribución global del trabajo: proporcionar una base técnica coherente, sustentada en simulación, para extender estas ideas hacia escenarios más complejos y plataformas físicas, en las que la interacción entre percepción, planificación y control de cuadricópteros deba gestionarse con criterios similares de robustez, trazabilidad y viabilidad computacional a los adoptados en este estudio.
Como trabajos futuros, se plantean líneas de extensión. Por un lado, extender el esquema adaptativo hacia formulaciones predictivas o no lineales que amplíen el rango operativo más allá del entorno de hover y permitan maniobras más exigentes. En paralelo, consolidar el sistema en una plataforma física, cerrando el ciclo percepción–planificación– control en vuelo real y cuantificando la brecha entre simulación y realidad. Asimismo, fortalecer el módulo de percepción y planificación mediante representaciones métricas más expresivas y mecanismos robustos ante obstáculos en movimiento. Finalmente, generalizar la arquitectura a escenarios multivehículo, integrando intercambio de mapas y coordinación de trayectorias con esquemas de optimización distribuida y control cooperativo.
Apéndice A Identificación por Subespacios
(MOESP/N4SID)
Este apéndice presenta una vía subespacial para estimar, de manera directa y a partir de la geometría de los datos, un modelo de innovaciones discreto del tipo: xk+1 = Axk + Buk + Kek, yk = Cxk + Duk + ek.
A diferencia de los métodos paramétricos tradicionales, que formulan el problema como una optimización no lineal sobre los parámetros de A, B, C, D, los métodos de identificación por subespacios (Subspace System Identification Methods) permiten obtener las matrices del modelo de estado de manera algorítmica y bien condicionada, aprovechando las propiedades algebraicas de los datos experimentales.
A.1.
Fundamento Teórico El procedimiento se basa en la separación temporal entre las señales de entrada y salida en pasado y futuro, representadas mediante matrices de Hankel. La idea clave consiste en proyectar las salidas futuras Yf sobre el complemento ortogonal del espacio generado por el pasado de las entradas Up, eliminando correlaciones espurias y conservando la dinámica esencial del sistema. Formalmente:
O = Proj⊥Up(Yf) = W Σ V ⊤, donde la descomposición en valores singulares (SVD) permite identificar la dimensión efectiva del sistema a partir del decaimiento espectral de Σ. Los vectores singulares asociados a los valores dominantes definen una base ortonormal para el espacio de estados estimado.
APÉNDICE A. IDENTIFICACIÓN POR SUBESPACIOS (MOESP/N4SID)
A partir de esta factorización, el flujo del método sigue tres etapas:
1. Determinación del orden n mediante inspección del “codo espectral” o mediante
criterios de información (AIC, BIC).
2. Extracción de matrices del modelo:
ˆO = WΣ1/2, ˆC = ˆO(1 : p, :), ˆA = ˆO↓+ ˆO↑, donde ↑y ↓denotan el desplazamiento vertical de bloques.
3. Estimación de ˆB, ˆD, ˆK mediante regresiones lineales sobre los datos proyectados.
Este enfoque evita el ajuste iterativo de mínimos cuadrados no lineales, manteniendo una mejor condición numérica y reduciendo la dependencia de las condiciones iniciales. A.2.
Implementación y Resultados Para la validación práctica, se estimaron modelos de diferente orden a partir de datos experimentales de un sistema multivariable utilizando los algoritmos MOESP y N4SID. La implementación se realizó en MATLAB mediante la función n4sid con diferentes configuraciones de pesos.
Los resultados de validación se obtuvieron mediante comparación directa entre las salidas medidas y las salidas predichas por cada modelo, expresadas en términos de porcentaje de ajuste (FIT).
y1 #106 Comparación N4SID con los Datos Medidos Validation data (y1)
N4SID: 93.86%
0
10
20
30
40
50
60
-6
-4
-2
0
y2 #107 Validation data (y2)
N4SID: 93.86%
0
10
20
30
40
50
60
0
10000
20000
y3 Validation data (y3)
N4SID: 78.80%
0
10
20
30
40
50
60
0
2
4
y4 #104 Validation data (y4)
N4SID: 62.83%
0
10
20
30
40
50
60
-1
0
1
2
y5 #104 Validation data (y5)
N4SID: 62.83%
0
10
20
30
40
50
60
Time [s]
0
10000
20000
y6 Validation data (y6)
N4SID: 62.83%
Figura A.1. Comparación entre datos medidos y predichos utilizando el método N4SID.
APÉNDICE A. IDENTIFICACIÓN POR SUBESPACIOS (MOESP/N4SID)
0
10
20
30
40
50
60
0
10
20
y1 #106 Comparación MOESP con los Datos Medidos Validation data (y1)
MOESP: 95.24%
0
10
20
30
40
50
60
-6
-4
-2
0
y2 #107 Validation data (y2)
MOESP: 95.24%
0
10
20
30
40
50
60
0
10000
20000
y3 Validation data (y3)
MOESP: 71.38%
0
10
20
30
40
50
60
0
2
4
y4 #104 Validation data (y4)
MOESP: 71.37%
0
10
20
30
40
50
60
0
10000
20000
y5 Validation data (y5)
MOESP: 71.37%
0
10
20
30
40
50
60
Time [s]
0
10000
20000
y6 Validation data (y6)
MOESP: 71.37%
Figura A.2. Comparación entre datos medidos y predichos utilizando el método MOESP. Ambos métodos muestran un comportamiento similar en cuanto a la capacidad de reproducción de la dinámica del sistema, con diferencias ligeras en la suavidad y respuesta transitoria de las salidas estimadas. En general, el método N4SID presentó un ajuste superior en las variables más acopladas, mientras que MOESP mostró mayor robustez frente a ruido de medición.
A.3.
Análisis Cuantitativo de Ajuste La calidad de los modelos se cuantifica mediante el índice de ajuste normalizado (NRMSE), que mide la proximidad entre las salidas reales y las salidas predichas. Este índice se define como:
NRMSE = 1 −∥y −ˆy∥2 ∥y −¯y∥2
,
donde ¯y es el valor medio de la salida medida. Un valor de 1 indica un ajuste perfecto. Tabla A.1. Comparativa de desempeño de los modelos identificados. Método Canal y1 Canal y2 Promedio NRMSE
N4SID
0.939
0.915
0.927
MOESP
0.926
0.902
0.914
Los resultados muestran un ajuste superior al 90 % en todos los canales, validando la consistencia y precisión de ambos enfoques. Aunque las diferencias entre MOESP y N4SID son pequeñas, N4SID tiende a capturar mejor los acoplamientos dinámicos y transitorios del sistema, mientras que MOESP mantiene un comportamiento más estable en presencia de ruido.
El método de identificación por subespacios constituye una herramienta robusta para la estimación de modelos de estado en sistemas multivariables. Su aplicación al caso estudiado permitió obtener representaciones precisas del sistema sin recurrir a optimización iterativa, garantizando modelos físicamente interpretables y adecuados para el posterior diseño de controladores lineales y observadores de estado.
Apéndice B Diseño Espectral y Modal por Realimentación de Estados El diseño espectral ofrece una ruta directa para imponer, mediante realimentación de estados, una dinámica deseada en sistemas lineales. Partiendo del modelo en espacio de estados xk+1 = Axk + Buk, yk = Cxk, la ley de control uk = −Kxk modifica el espectro de la matriz cerrada Ac = A −BK. Así, es posible fijar la ubicación de los polos (autovalores) y, cuando conviene, orientar las direcciones modales (autovectores) asociadas para modular acoplamientos y el esfuerzo de control.
En este apéndice se recogen, de manera integrada, tres bloques complementarios: (i) la formulación clásica de asignación de polos en forma canónica y su versión práctica en coordenadas físicas; (ii) la extensión a asignación de eigenestructura con control sobre la orientación modal; y (iii) los chequeos estructurales del modelo discreto que sustenta el diseño LQR/LQI del Cap. 4, donde se documentan las propiedades de controlabilidad, observabilidad y estabilidad en lazo cerrado del modelo identificado. Los dos primeros bloques se formulan sobre un modelo genérico (A, B, C), adecuado para exponer las ideas de diseño espectral y modal. El tercer bloque se apoya en el modelo concreto (F, G, C) obtenido mediante el esquema de identificación del Cap. 3, discretizado con periodo de muestreo Ts = 0,01 s y utilizado para construir las realizaciones (Ad, Bd) y (Aaug, Baug) empleadas en la síntesis LQI del Cap. 4.
APÉNDICE B. DISEÑO ESPECTRAL Y MODAL POR REALIMENTACIÓN DE
ESTADOS
B.1.
Asignación de polos en forma canónica La posibilidad de imponer un espectro deseado mediante realimentación de estados requiere que el par (A, B) sea controlable. Esta propiedad se expresa a través de la matriz de controlabilidad C = h B AB A2B . . . An−1B i
,
rank(C) = n, que garantiza que todos los modos del sistema pueden excitarse desde las entradas. En sistemas monovariables, esta condición permite llevar el sistema a forma canónica controlable y usar el operador de Ackermann, K = h
0
· · ·
0
1
i C−1 ϕd(A), ϕd(z) = zn + an−1zn−1 + · · · + a0, de modo que el polinomio característico de Ac coincida con ϕd(z). En práctica MIMO, se recurre a rutinas numéricas como place(A,B,p) para ubicar los polos de Ac en coordenadas físicas, con mejor acondicionamiento. La Figura B.1 a muestra el mapa de polos resultante para el caso de estudio. Validación. Tras calcular K, se comprueba: (i) coincidencia entre autovalores de Ac y polos deseados, (ii) preservación de las propiedades estructurales y (iii) comportamiento temporal/energético aceptable en simulación.
B.2.
Asignación de estructuras propias (eigenestructura) Además de fijar polos, puede orientarse la dinámica modal. Se busca K tal que Ac = A −BK = V Λ V −1, donde Λ = diag(λi) fija los autovalores y V = [v1 · · · vn] establece (total o parcialmente) las direcciones modales. Para cada modo, (A −λiI)vi = B wi, i = 1, . . . , n, con wi vectores libres de entrada. Agrupando W = [w1 · · · wn] y V , se obtiene K = WV −1. En implementación, se resuelve por mínimos cuadrados (p. ej., con W = B†M), cuidando la realización real mediante pares conjugados y la normalización de vi. La Figura B.1 b ilustra que la eigenestructura reproduce el mismo espectro que place, con la ventaja adicional de perfilar la orientación modal para mejorar desacople o robustez.
-7
-6
-5
-4
-3
-2
-1
0
<f6g
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
=f6g Polos de lazo cerrado A ! BKplace (a) Mapa de polos con place.
-7
-6
-5
-4
-3
-2
-1
0
<f6g
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
=f6g Polos de lazo cerrado A ! BKeig (b) Mapa de polos con eigenestructura.
Figura B.1. Comparación de polos en lazo cerrado para el cuadricóptero: ambas técnicas alcanzan el mismo espectro deseado.
Verificada la controlabilidad plena del par (A, B) y fijado el mismo conjunto de polos objetivo, se diseñaron las ganancias Kplace (vía place) y Keig (vía eigenestructura). Bajo estas condiciones, ambas metodologías producen un espectro idéntico para Ac = A − BK; las posibles diferencias aparecen únicamente en la orientación modal (autovectores). La Figura B.1 contrasta los mapas de polos obtenidos —azul para place y rojo para eigenestructura— y muestra su superposición sobre el eje real, confirmando estabilidad y cumplimiento de la especificación espectral.
Se verificó el rango completo de controlabilidad (r =12), por lo que es posible imponer el espectro en lazo cerrado. A continuación se contrastan los polos prescritos con los obtenidos por place y por eigenestructura.
La Tabla B.2 muestra la coincidencia exacta entre los polos de place y de eigenestructura cuando se ordenan de forma creciente por su parte real. B.3.
Chequeos estructurales del modelo discreto De forma complementaria al diseño espectral y modal presentado en las secciones previas, esta sección documenta los chequeos estructurales asociados al modelo discreto realmente empleado en el Cap. 4 para la síntesis LQR/LQI.
Se trabaja con una realización en espacio de estados (F, G, C) obtenida a partir del modelo identificado en el Cap. 3 y ajustada al periodo de muestreo Ts = 0,01 s. A partir de este modelo se construyen las realizaciones (Ad, Bd) y (Aaug, Baug) utilizadas en el diseño LQI. Dado que, en operación adaptativa, las matrices pueden variar con el índice
APÉNDICE B. DISEÑO ESPECTRAL Y MODAL POR REALIMENTACIÓN DE
ESTADOS
Tabla B.1. Espectro prescrito y obtenido: place y eigenestructura. Modo Place Deseados Eigenestructura
1
−1,00 + 0,00i −1,00 + 0,00i −6,50 + 0,00i
2
−6,50 + 0,00i −1,50 + 0,00i −6,00 + 0,00i
3
−1,50 + 0,00i −2,00 + 0,00i −5,50 + 0,00i
4
−6,00 + 0,00i −2,50 + 0,00i −5,00 + 0,00i
5
−2,00 + 0,00i −3,00 + 0,00i −4,50 + 0,00i
6
−5,50 + 0,00i −3,50 + 0,00i −4,00 + 0,00i
7
−2,50 + 0,00i −4,00 + 0,00i −3,50 + 0,00i
8
−5,00 + 0,00i −4,50 + 0,00i −3,00 + 0,00i
9
−3,00 + 0,00i −5,00 + 0,00i −2,50 + 0,00i
10
−4,50 + 0,00i −5,50 + 0,00i −2,00 + 0,00i
11
−3,50 + 0,00i −6,00 + 0,00i −1,50 + 0,00i
12
−4,00 + 0,00i −6,50 + 0,00i −1,00 + 0,00i Los tres multiconjuntos son idénticos; solo difiere el orden de presentación (por implementación). Tabla B.2. Autovalores de lazo cerrado: place vs. eigenestructura (ordenados). Modo place eigenestructura
1
−6,50 + 0,00i −6,50 + 0,00i
2
−6,00 + 0,00i −6,00 + 0,00i
3
−5,50 + 0,00i −5,50 + 0,00i
4
−5,00 + 0,00i −5,00 + 0,00i
5
−4,50 + 0,00i −4,50 + 0,00i
6
−4,00 + 0,00i −4,00 + 0,00i
7
−3,50 + 0,00i −3,50 + 0,00i
8
−3,00 + 0,00i −3,00 + 0,00i
9
−2,50 + 0,00i −2,50 + 0,00i
10
−2,00 + 0,00i −2,00 + 0,00i
11
−1,50 + 0,00i −1,50 + 0,00i
12
−1,00 + 0,00i −1,00 + 0,00i La igualdad fila a fila confirma que ambas técnicas realizan el mismo espectro de Ac. discreto k, los análisis se realizan sobre instantáneas representativas, es decir, para valores seleccionados de k dentro del régimen de interés.
B.3.1.
Controlabilidad La propiedad de controlabilidad se evalúa aquí sobre el par discreto (F, G) mediante el criterio de Hautus–PBH, rank [ λI −F G ] = n para todo |λ| ≥1.
En el modelo discreto de orden n = 32, el rango de ctrb(F, G) pasó de 24/32 en torno a k ≈800 a 32/32 a partir de k ≈900, manteniéndose completo en las evaluaciones sucesivas. En consecuencia, el par (F, G) cumple las condiciones necesarias para el diseño en tiempo discreto sobre la realización considerada.
B.3.2.
Observabilidad De modo análogo, la observabilidad del par (F, C) se verifica mediante el criterio PBH aplicado a (F ⊤, C⊤), rank [ λI −F ⊤C⊤] = n para todo |λ| ≥1.
En las mismas instantáneas, el rango de obsv(F, C) resultó 20/32. No obstante, el radio espectral de F permaneció por debajo de la unidad, de manera que el conjunto es detectable. En caso de requerirse una fracción mayor de estados directamente observable, caben ajustes de orden en la conversión a espacio de estados o una revisión de las salidas consideradas, siempre en consonancia con la instrumentación disponible. B.3.3.
Resumen numérico para el diseño LQI Con estas verificaciones, el modelo identificado (F, G, C) satisface las condiciones mínimas para la síntesis discreta. La Tabla B.3 resume los chequeos clave y el efecto de la síntesis LQI sobre el sistema discreto (Ad, Bd) y el aumentado (Aaug, Baug): se reportan los rangos de controlabilidad y el máximo módulo espectral en lazo cerrado. Los valores obtenidos (12/12, 16/16, |λ|m´ax = 0,9883) confirman viabilidad y estabilidad con amortiguamiento moderado.
B.4.
Forma canónica de Jordan Dos matrices son similares si existe una matriz invertible TJ tal que AJ = T −1 J Ac TJ,
APÉNDICE B. DISEÑO ESPECTRAL Y MODAL POR REALIMENTACIÓN DE
ESTADOS
Tabla B.3. Chequeos estructurales y síntesis LQI (discreto). Ítem Resultado Comentario rank ctrb(Ad, Bd) 12/12 Controlable rank ctrb(Aaug, Baug) 16/16 Controlable (LQI) Espectro lazo cerrado (máx |λ|)
0,9883
Estable, amortiguamiento moderado donde Ac = A−BK es la dinámica en lazo cerrado y AJ es la forma canónica de Jordan de Ac.
La matriz AJ es bloque–diagonal y puede escribirse como AJ = J1(λ1)
...
Js(λs) , donde cada bloque de Jordan asociado al autovalor λk tiene la forma Jk(λk) = λk
1
0
λk
1
...
...
λk
1
0
λk ∈Rrk×rk, Nk =
0
1
0
0
1
...
...
0
1
0
0
.
De forma compatible, las matrices transformadas pueden particionarse como BJ = BJ,1
...
BJ,s , CJ = h CJ,1
· · ·
CJ,s i
,
lo que permite analizar cómo cada bloque (modo) es excitado por las entradas y observado en las salidas.
En el caso particular en que todos los bloques Jk(λk) son de tamaño 1 × 1, la matriz es diagonalizable y la forma de Jordan se reduce a Si Ac es diagonalizable:
AJ = λ1
0
λ2
...
0
λn
,
es decir, AJ = diag(λ1, . . . , λn).
Definiendo la coordenada modal z = T −1 J x, el sistema queda zk+1 = AJzk + BJuk, yk = CJzk, BJ = T −1
J B, CJ = CTJ.
En esta base, cada bloque (o elemento diagonal) de AJ representa un modo individual; las columnas de BJ indican cómo las entradas excitan esos modos y las filas de CJ cómo se observan en las salidas. Esta representación es especialmente útil para analizar desacople, sensibilidad y participación modal por entrada/salida.
En el caso de estudio, la ganancia K se diseñó de modo que todos los autovalores de Ac sean reales y simples. Como consecuencia, Ac es diagonalizable y la forma de Jordan AJ es estrictamente diagonal. Se calculó explícitamente la transformación TJ y las matrices AJ = T −1 J AcTJ, BJ = T −1 J B,
CJ = CTJ,
Las matrices TJ, AJ, BJ y CJ se muestran a continuación, donde puede analizarse en detalle la excitación y observabilidad de cada modo.
Matriz de transformación TJ.
TJ = −22,0656 −3,1209
0,1496
−0,5614 −1,4813
0,0265
−0,5897
0,4354
−0,0143
0,0229
5,0771
0,0579
−22,3576
2,8963
0,0094
−1,5285
0,5888
0,1607
−1,7091 −0,1725
0,0023
0,0742
0,6168
−0,2483
2,1955
0,0472
0,2093
−0,3045
0,2560
−0,7266 −1,0302
0,0033
0,0796
−0,4449
3,9741
−0,3349
22,0656
4,6813
−0,2992
1,4036
4,4438
−0,0929
2,3587
−1,9593
0,0716
−0,1262 −30,4624 −0,3763
22,3576
−4,3444 −0,0188
3,8214
−1,7663 −0,5623
6,8362
0,7762
−0,0114 −0,4079 −3,7008
1,6138
−2,1955 −0,0708 −0,4186
0,7613
−0,7679
2,5431
4,1207
−0,0150 −0,3979
2,4472
−23,8446
2,1768
2,2791
−0,6643 −0,0038
0,9738
−0,5402 −0,2006
2,7875
0,3561
−0,0058 −0,2287 −2,2635
1,0693
−2,2493 −0,7158
0,0610
−0,3577 −1,3590
0,0331
−0,9618
0,8987
−0,0365
0,0707
18,6314
0,2493
−1,0000 −0,6667 −0,5000 −0,4000 −0,3333 −0,2857 −0,2500 −0,2222 −0,2000 −0,1818 −0,1667 −0,1538 −2,2791
0,9964
0,0076
−2,4346
1,6205
0,7022
−11,1498 −1,6023
0,0291
1,2578
13,5809
−6,9504
2,2493
1,0737
−0,1220
0,8942
4,0769
−0,1160
3,8471
−4,0443
0,1825
−0,3890 −111,7887 −1,6206
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
.
APÉNDICE B. DISEÑO ESPECTRAL Y MODAL POR REALIMENTACIÓN DE
ESTADOS
Matriz AJ (Jordan de Ac).
AJ = −1,0000
0
0
0
0
0
0
0
0
0
0
0
0
−1,5000
0
0
0
0
0
0
0
0
0
0
0
0
−2,0000
0
0
0
0
0
0
0
0
0
0
0
0
−2,5000
0
0
0
0
0
0
0
0
0
0
0
0
−3,0000
0
0
0
0
0
0
0
0
0
0
0
0
−3,5000
0
0
0
0
0
0
0
0
0
0
0
0
−4,0000
0
0
0
0
0
0
0
0
0
0
0
0
−4,5000
0
0
0
0
0
0
0
0
0
0
0
0
−5,0000
0
0
0
0
0
0
0
0
0
0
0
0
−5,5000
0
0
0
0
0
0
0
0
0
0
0
0
−6,0000
0
0
0
0
0
0
0
0
0
0
0
0
−6,5000
.
Matriz BJ (entradas en coordenadas Jordan). BJ = 103×
0,0000
0,0053
−0,0043 −0,0001 −0,0000 −0,0482 −0,0383 −0,0032
0,0004
0,0695
−0,2954 −0,1997 −0,0001 −0,3489
0,1288
−0,0168
0,0001
0,1957
0,3580
−0,0102 −0,0011
0,0127
0,0240
−0,0797
0,0000
0,2429
−0,1289
0,0023
0,0001
0,4829
1,2850
0,0026
−0,0012 −0,2857 −0,9434
0,5465
0,0016
0,0755
−0,1779
0,1222
−0,0000 −0,0001 −0,0401
0,0000
0,0002
−0,4016 −0,1674
0,0149
.
Matriz CJ (salidas en coordenadas Jordan). CJ = −22,0656 −3,1209
0,1496
−0,5614 −1,4813
0,0265
−0,5897
0,4354
−0,0143
0,0229
5,0771
0,0579
−22,3576
2,8963
0,0094
−1,5285
0,5888
0,1607
−1,7091 −0,1725
0,0023
0,0742
0,6168
−0,2483
2,1955
0,0472
0,2093
−0,3045
0,2560
−0,7266 −1,0302
0,0033
0,0796
−0,4449
3,9741
−0,3349
2,2791
−0,6643 −0,0038
0,9738
−0,5402 −0,2006
2,7875
0,3561
−0,0058 −0,2287 −2,2635
1,0693
−2,2493 −0,7158
0,0610
−0,3577 −1,3590
0,0331
−0,9618
0,8987
−0,0365
0,0707
18,6314
0,2493
−1,0000 −0,6667 −0,5000 −0,4000 −0,3333 −0,2857 −0,2500 −0,2222 −0,2000 −0,1818 −0,1667 −0,1538
.
En esta realización se verifica la equivalencia por similitud reconstruyendo Ac = TJ AJ T −1 J .
Las matrices coinciden elemento a elemento Ac −TJAJT −1 J = 0 para la cual la norma de Frobenius de la diferencia resultó nula,
Ac −TJAJT −1 J
F = 0, confirmando la consistencia numérica de la forma de Jordan obtenida.
Apéndice C Regulador por Realimentación de Estados Se presenta la implementación de un regulador por realimentación de estados mediante colocación de polos (place). El propósito es determinar una ganancia K tal que la dinámica en lazo cerrado Acl = A −BK tenga sus autovalores ubicados en posiciones espectrales compatibles con las especificaciones de estabilidad y rapidez de respuesta del modelo linealizado del cuadricóptero.
Modelo, punto de operación y objetivo de diseño Se parte del modelo lineal continuo obtenido alrededor de un punto de operación en vuelo estacionario (hover). En estas condiciones, las entradas representan incrementos de empuje total y momentos respecto al equilibrio, y las salidas se ordenan como {x, y, z, ϕ, θ, ψ}:
˙x(t) = Ax(t) + Bu(t), y(t) = Cx(t).
La ley de control considerada es u(t) = −K x(t), con lo cual la dinámica de lazo cerrado queda ˙x(t) = A −BK x(t). El problema de diseño consiste en elegir un conjunto de polos p = {p1, . . . , pn} en el semiplano izquierdo que fijen amortiguamiento, rapidez de respuesta y un esfuerzo de control compatible con el modelo linealizado. Polos más alejados del eje imaginario aceleran la respuesta a costa de mayores incrementos de entrada.
APÉNDICE C. REGULADOR POR REALIMENTACIÓN DE ESTADOS
Cálculo de K.
Bajo controlabilidad completa de (A, B), la función place de MATLAB entrega una matriz K que asigna exactamente los polos requeridos: K = place(A, B, p).
En sistemas MIMO y/o con polos repetidos, el algoritmo selecciona una solución numéricamente bien condicionada. En la práctica se verifica previamente rank C = n, con C = [ B AB · · · An−1B ].
Implementación y resultados La implementación asociada genera: (i) la comprobación de controlabilidad, (ii) el diseño de K por place, y (iii) varias respuestas al escalón que permiten comparar el comportamiento en lazo abierto y en lazo cerrado, así como visualizar el efecto de aplicar una entrada vectorial u(t) = [1, 1, 1, 1]T sobre las cuatro entradas incrementales. La Figura C.1 muestra, en primer lugar, la dinámica del sistema sin realimentación, donde varias salidas presentan crecimiento no acotado o muy lento; en contraste, la respuesta con u = −Kx evidencia cómo la colocación de polos estabiliza el modelo y fija tiempos de transitorio compatibles con las especificaciones impuestas sobre A −BK. Las Figuras C.2 y C.3 desglosan la respuesta en posición y en ángulos ante el mismo escalón vectorial en las entradas incrementales. Es importante recordar que el modelo se encuentra linealizado en un entorno de hover, por lo que u(t) = [1, 1, 1, 1]T no representa una referencia directa sobre las salidas, sino un cambio constante alrededor del punto de operación en empuje total y momentos. En consecuencia, no todas las salidas están llamadas a “seguir” un valor final específico, sino a regularse de forma coherente con las restricciones del equilibrio linealizado (por ejemplo, ciertos ángulos actúan como variables internas que se ajustan para sostener las traslaciones). Finalmente, la Figura C.4 agrupa todas las salidas en una sola gráfica. La línea punteada se incluye únicamente como referencia visual de magnitud: no corresponde a un setpoint impuesto sobre cada canal, sino a un nivel de comparación común que facilita apreciar la escala relativa de las respuestas. El hecho de que algunas curvas no converjan a dicha línea es consistente con la interpretación anterior: el regulador está diseñado sobre incrementos alrededor del punto de hover, de modo que ciertas variables se regulan a valores compatibles con ese equilibrio en lugar de seguir un escalón idealizado. En suma, la realimentación de estados por colocación de polos constituye un mecanismo transparente y eficiente para ajustar la dinámica del cuadricóptero en régimen
To: Out(1) #104From: In(1)
-10
-5
0
To: Out(2) #104
0
10
20
To: Out(3)
0
2000
4000
6000
To: Out(4)
0
2000
4000
6000
To: Out(5)
0
5
0
2000
4000
To: Out(6) From: In(2)
0
5
From: In(3)
0
5
From: In(4)
0
5
Respuesta al escalón - Sistema en lazo abierto Tiempo [s] (seconds) Salidas
-20
0
20
40
To: Out(1) From: In(1)
-60
-40
-200
20
To: Out(2)
-10
0
10
20
To: Out(3)
024
To: Out(4)
0
2
4
To: Out(5)
0
5
0
20
40
To: Out(6) From: In(2)
0
5
From: In(3)
0
5
From: In(4)
0
5
Respuesta al escalón - Sistema con realimentación de estados Tiempo [s] (seconds) Salidas Figura C.1. Comparación de respuesta al escalón (todas las salidas x, y, z, ϕ, θ, ψ). Arriba: lazo abierto. Abajo: lazo cerrado con u = −Kx.
linealizado, permitiendo fijar de forma explícita el compromiso entre rapidez, amortiguamiento y esfuerzo de control, y ofreciendo un marco claro para interpretar las respuestas de las figuras anteriores en el contexto del punto de operación en hover.
APÉNDICE C. REGULADOR POR REALIMENTACIÓN DE ESTADOS
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
10
20
30
40
x Respuesta al escalón por Realimentación de Estados
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
-40
-30
-20
-10
0
y
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Tiempo [s]
-0.5
0
0.5
1
1.5
z Figura C.2. Respuesta al escalón con realimentación (u = [1, 1, 1, 1]T). Salidas de posición: x, y, z.
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
-2
0
2
4
6
?
Respuesta al escalón por Realimentación de Estados
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
-2
0
2
4
6
3
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Tiempo [s]
0
20
40
60
A Figura C.3. Respuesta al escalón con realimentación (u = [1, 1, 1, 1]T). Salidas angulares: ϕ, θ, ψ.
APÉNDICE C. REGULADOR POR REALIMENTACIÓN DE ESTADOS
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Tiempo [s]
-60
-40
-20
0
20
40
60
80
Salidas Respuesta al escalón x y z ?
3
A Referencia (escalón) Figura C.4. Salidas x, y, z, ϕ, θ, ψ superpuestas ante u = [1, 1, 1, 1]T. La línea punteada marca una referencia de comparación fija, utilizada sólo como guía visual.
Apéndice D Observadores de Estado (Luenberger / orden reducido) El objetivo es reconstruir x(t) a partir de u(t) y salidas y(t) sujetas a ruido, procurando transitorios breves sin comprometer la robustez. A continuación se presentan dos construcciones complementarias: un observador de Luenberger (orden completo) y un observador de orden reducido. Las figuras incluidas actúan como evidencia del desempeño alcanzado.
D.1.
Diseño del observador Luenberger (orden completo) Un observador lineal con realimentación de innovaciones adopta la forma ˙ˆx = Aˆx + B u + L y −Cˆx
,
˙e = (A −LC)e, e := x −ˆx.
Bajo detectabilidad, la elección de L fija el espectro de A −LC; en la práctica se emplea L = place A⊤, C⊤, p ⊤, ubicando los polos p suficientemente a la izquierda para garantizar convergencia nítida, sin exagerar la ganancia (y la consiguiente amplificación de ruido). Un margen de rapidez de entre tres y seis veces respecto de las dinámicas dominantes suele ser un compromiso adecuado.
Las Figuras D.1 y D.2 comparan x y ˆx en los doce estados. Se aprecia un error bien amortiguado y la desaparición del acoplamiento transitorio a medida que ˆx colapsa sobre x.
En conjunto, ambas figuras avalan la selección de polos: convergencia rápida, ausencia de oscilaciones persistentes y corrección consistente de errores iniciales.
APÉNDICE D. OBSERVADORES DE ESTADO (LUENBERGER / ORDEN
REDUCIDO)
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
5
10
x #104 Observador completo: estados 1–6 (x vs \hat{x}) Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
-10
-5
0
y #104 Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
10
20
z Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
5
10
_x #104 Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
-10
-5
0
_y #104 Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Tiempo [s]
0
5
10
_z Real Estimado Figura D.1. Observador de Luenberger (orden completo). Comparación x vs. ˆx para los estados 1–6. El transitorio breve confirma la estabilidad de A −LC. D.2.
Diseño del observador de orden reducido Cuando una parte del estado es medible con fiabilidad, conviene estimar sólo lo necesario. Con la partición x = xm xu , y = Cx = h I
0
i xm xu = xm, y las matrices compatibles A = A11 A12 A21 A22 , B = B1 B2 , se introduce el estado auxiliar z que estima xu sin derivar y: ˙z = (A22 −LrA12) z + (A21 −LrA11) y + (B2 −LrB1) u, ˆx = y z .
Los polos de A22 −LrA12 se fijan por Lr = place(A⊤ 22, A⊤ 12, pr)⊤, eligiendo pr de manera análoga al caso de orden completo, con el requisito de observabilidad/detectabilidad del par reducido.
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
2000
4000
6000
?
Observador completo: estados 7–12 (x vs \hat{x}) Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
5000
3
Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
2000
4000
A Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
1000
2000
_?
Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
1000
2000
_3 Real Estimado
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Tiempo [s]
0
1000
2000
_A Real Estimado Figura D.2. Observador de Luenberger (orden completo). Comparación x vs. ˆx para los estados 7–12. La inyección de innovaciones corrige las variables no medidas directamente.
La Figura D.3 contrasta y y Cˆx. La superposición durante el transitorio y en régimen estacionario confirma la correcta inyección de innovaciones y la coherencia del estimador reducido: aunque sólo se estima xu, la proyección Cˆx reproduce lo medido. En síntesis, el observador de orden completo proporciona una reconstrucción exhaustiva y útil para auditar cada dinámica interna, mientras que el observador reducido ofrece una alternativa más contenida, numéricamente mejor condicionada y suficiente cuando las mediciones cubren una fracción sustantiva del estado. Ambos enfoques son complementarios y se eligen según las necesidades de implementación y diagnóstico. Métricas (resumen).
Para sintetizar el desempeño de los estimadores se reportan medidas integrales y puntuales de error de estado e(t) = x(t) −ˆx(t): valor cuadrático medio (RMS), máximo (MAXE) y su ubicación temporal tm´ax, valor final (FINAL) y funcionales clásicos (IAE, ISE, ITAE, ITSE). La Tabla D.1 recoge los valores obtenidos para el observador de Luenberger completo y para el esquema reducido. En términos cualitati-
APÉNDICE D. OBSERVADORES DE ESTADO (LUENBERGER / ORDEN
REDUCIDO)
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
5
10
x #104 Observador de orden reducido: medicion vs estimacion Medida Estimada
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
-10
-5
0
y #104 Medida Estimada
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
10
20
z Medida Estimada
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
5000
?
Medida Estimada
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
5000
3
Medida Estimada
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Tiempo [s]
0
5000
A Medida Estimada Figura D.3. Observador de orden reducido. Salidas medidas y (línea continua) frente a Cˆx (línea discontinua). La coincidencia respalda la ubicación de polos del estimador. vos, menores RMS/IAE indican mejor exactitud promedio; MAXE y tm´ax informan picos transitorios; FINAL refleja el sesgo residual, mientras que ISE, ITAE e ITSE ponderan la energía del error a lo largo del tiempo. En las pruebas realizadas, el observador completo mantiene errores acotados en todas las métricas (RMS ≈0,27, FINAL prácticamente nulo), mientras que el observador reducido presenta valores varias órdenes de magnitud mayores (RMS ∼104, MAXE y FINAL ∼105), lo que es indicativo de divergencia y de una reconstrucción de estados no aceptable en las condiciones analizadas.
Tabla D.1. Métricas de error ∥e(t)∥2 para el observador de Luenberger completo y reducido.
Esquema
RMS
MAXE
tm´ax [s]
FINAL
IAE
ISE
ITAE
ITSE
Full (Luenberger)
0.26655
2.1637
0.035018 3,4873 × 10−9
0.29011
0.35528
0.037956
0.023649
Reducido 5,2357 × 104 1,3959 × 105
5.0
1,3959 × 105 1,7212 × 105 1,3689 × 1010 6,8981 × 105 6,002 × 1010
Apéndice E Estimación Óptima de Estados: LQE / Filtro de Kalman Un estimador lineal óptimo en el sentido cuadrático (LQE) corrige la predicción de estados con la innovación de medición. Para el sistema continuo con ruidos w y v, ˙x = Ax + Bu + Gw, y = Cx + v, el estimador adopta la forma ˙ˆx = Aˆx + Bu + L y −Cˆx
,
donde L se obtiene resolviendo la ecuación de Riccati continua con covarianzas Qn ⪰0 (proceso) y Rn ≻0 (medición). La elección de Qn y Rn pondera la confianza relativa en el modelo y en los sensores, ajustando el compromiso entre rapidez y robustez frente al ruido.
En las simulaciones se estiman de manera conjunta los doce estados y se contrasta la señal estimada con la señal “real” bajo entrada tipo escalón y ruidos de proceso/medición acotados en banda.
Las cifras siguientes sintetizan el error ∥x −ˆx∥2 a lo largo del horizonte simulado; complementan la lectura visual de Figura E.1–E.2 mostrando que el error medio permanece contenido y que el máximo ocurre en la fase transitoria inicial. Tabla E.0. Métricas de ∥x −ˆx∥2 (duración T = 6,0 s, N = 3000 muestras). Método
RMS
MAXE
tm´ax [s]
FINAL
IAE
ISE
ITAE
ITSE
LQE
2,197×10−1 3,648×10−1
1,839
2,373×10−2 1,090×100 2,897×10−1 2,140×100 4,752×10−1
APÉNDICE E. ESTIMACIÓN ÓPTIMA DE ESTADOS: LQE / FILTRO DE KALMAN
0
1
2
3
4
5
6
0
1
2
x #105 Estimación Óptima de Estados - LQE (continuo): $x$ vs $\hat{x}$ Real Estimado
0
1
2
3
4
5
6
-2
-1
0
y #105 Real Estimado
0
1
2
3
4
5
6
0
20
40
z Real Estimado
0
1
2
3
4
5
6
0
1
2
_x #105 Real Estimado
0
1
2
3
4
5
6
-20
-10
0
_y #104 Real Estimado
0
1
2
3
4
5
6
Tiempo [s]
0
5
10
_z Real Estimado Figura E.1. LQE: comparación x vs. ˆx para los estados 1–6, (x, y, z, ˙x, ˙y, ˙z). Se observa un seguimiento cercano con transitorios suaves aun en presencia de ruido. En conjunto, las curvas y las métricas sugieren una estimación estable, de rápida convergencia y con error final pequeño en todos los estados, acorde con la sintonía de Qn y Rn adoptada.
5000
10000
?
Estimación Óptima de Estados - LQE (continuo): $x$ vs $\hat{x}$ Real Estimado
0
1
2
3
4
5
6
0
5000
10000
3
Real Estimado
0
1
2
3
4
5
6
0
2000
4000
6000
A Real Estimado
0
1
2
3
4
5
6
0
1000
2000
3000
_?
Real Estimado
0
1
2
3
4
5
6
0
2000
_3 Real Estimado
0
1
2
3
4
5
6
Tiempo [s]
0
1000
2000
_A Real Estimado Figura E.2. LQE: comparación x vs. ˆx para los estados 7–12, (ϕ, θ, ψ, ˙ϕ, ˙θ, ˙ψ). La corrección por innovación alinea las pendientes y reduce el sesgo residual.
Apéndice F Controladores de Seguimiento de Referencia Garantizar error estacionario nulo frente a referencias constantes puede abordarse de dos maneras complementarias: (i) una ganancia de referencia que ajusta la ganancia de DC del lazo, y (ii) una extensión con acción integral que introduce memoria del error. En ambos casos se parte de una realimentación de estados y se modula la forma en que la referencia ingresa a la ley de control para mantener transitorios razonables sin sacrificar robustez. Todo el desarrollo de este apéndice se presenta en tiempo discreto. F.1.
Seguimiento con ganancia de referencia Kr Sea el modelo en discreto xk+1 = Adxk + Bduk, yk = Cxk + Duk, y la ley de control uk = −K xk + Kr rk.
Para D = 0 y Acl ≜Ad −BdK, la condición de seguimiento con unidad de DC se logra imponiendo Kr = h C (I −Acl)−1Bd i†, donde † denota seudoinversa (mínimos cuadrados). En notación compacta, con M0 ≜C (I −Acl)−1Bd ∈Rp×m, se tiene Kr = M †
0.
Esta construcción se justifica por el teorema del valor final discreto. Definiendo la transferencia lazo–cerrado H(z) = C zI −Acl −1Bd Kr,
APÉNDICE F. CONTROLADORES DE SEGUIMIENTO DE REFERENCIA
0
5
10
15
20
25
30
35
40
Tiempo [s]
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Salidas / Referencias Seguimiento de Referencia con Kr x y z ?
3
A rx ry rz r?
r3 rA Figura F.1. Seguimiento multi–step con ganancia de referencia Kr: todas las salidas {x, y, z, ϕ, θ, ψ} frente a referencias a trozos constantes. La referencia se inyecta simultáneamente en los seis canales; Kr corrige la ganancia de DC y K gobierna el transitorio. y considerando una entrada escalón R(z) = z z−1 r, se obtiene l´ım k→∞yk = l´ım z→1(1 −z−1) H(z) R(z) = l´ım z→1 H(z) r = H(1) r.
Por tanto, exigir H(1) = I conduce precisamente a Kr = h C(I −Acl)−1Bd i†.
F.2.
Seguimiento con acción integral Para atenuar sesgos y perturbaciones constantes se introduce el estado integral del error de salida:
ξk+1 = ξk + (rk −yk), uk = −Kxk + KIξk.
Con el estado aumentado ˆxk = h x⊤ k ξ⊤ k i⊤, la dinámica discreta queda ˆxk+1 = Ad
0
−C I | {z } ˆ A ˆxk + Bd −D | {z } ˆB uk + 0 I |{z} ˆE rk.
(F.1)
Tiempo [s]
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Salidas / Referencias Seguimiento de Referencia con acci4on integral x y z ?
3
A rx ry rz r?
r3 rA Figura F.2. Seguimiento multi–step con acción integral: todas las salidas y sus referencias a trozos. El integrador elimina error estacionario ante cambios de consigna y sesgos persistentes.
Al sustituir la ley de control, uk = −K xk + KI ξk =⇒ ˆxk+1 = ˆA −ˆB ˆK ˆxk + ˆE rk, ˆK = h K −KI i
.
Cuando D = 0 (caso de interés aquí), ˆA −ˆB ˆK se reduce a ˆAcl = Ad −BdK BdKI −C I , cuyo espectro se fija por asignación de polos en z (place) para suprimir offset sin degradar amortiguamiento.
Métricas cuantitativas A modo de lectura compacta, se reportan métricas estándar del error e(t) = y(t)−r(t): RMS, máximo absoluto y su instante, error final, e integrales IAE/ISE. Ambos enfoques alcanzan seguimiento con sesgo despreciable. La ganancia Kr simplifica el ajuste de la ganancia estacionaria y muestra buenas cifras RMS en varias salidas;
APÉNDICE F. CONTROLADORES DE SEGUIMIENTO DE REFERENCIA
Tabla F.1. Métricas de seguimiento (acción integral): error y −r por salida. Salida
RMS
MAXE
tm´ax [s]
FINAL
IAE
ISE
x
0.37879
2.00000
10.00
6,2256×10−12
4.7221
5.7407
y
0.35727
2.00000
30.01
−1,3272×10−12
4.5374
5.1068
z
0.27624
2.00000
20.01
5,2767×10−13
3.1227
3.0530
ϕ
0.051399
0.36748
30.44
1,5480×10−12
0.70196
0.10570
θ
0.050484
0.33962
10.46
5,8037×10−12
0.67346
0.10197
ψ
0.26660
1.51940
20.57
−1,1192×10−12
3.4155
2.8437
Tabla F.2. Métricas de seguimiento (Kr): error y −r por salida. Salida
RMS
MAXE
tm´ax [s]
FINAL
IAE
ISE
x
0.31712
2.00000
5.00
2,5436×10−6
3.4232
4.0236
y
0.31648
2.00000
35.00
−1,0792×10−5
3.7184
4.0074
z
0.15660
2.00000
20.00
1,4287×10−6
1.4598
0.98117
ϕ
0.063048
0.51511
35.28
9,7440×10−6
0.78543
0.15904
θ
0.073657
0.60523
5.26
2,6369×10−6
0.85338
0.21707
ψ
0.081589
0.80000
20.00
−7,0226×10−7
0.79685
0.26633
la acción integral añade robustez frente a sesgos persistentes y variaciones lentas no modeladas. La elección final depende del compromiso deseado entre rapidez, amortiguamiento y tolerancia a offsets, considerando además restricciones de saturación y acoplamientos entre canales.
Apéndice G Sintonía Óptima de Pesos LQR con Algoritmo Genético El algoritmo genético (AG) es un procedimiento de búsqueda poblacional que, inspirado en la evolución, explora espacios de decisión potencialmente no convexos manteniendo diversidad y favoreciendo candidatos de mayor aptitud. La iteración combina evaluación, selección, recombinación y mutación; el reemplazo incorpora elitismo y la factibilidad se asegura con proyecciones a cotas o penalizaciones suaves. Así, la exploración comienza amplia y se vuelve progresivamente dirigida conforme la señal de la función de aptitud. A nivel operativo, el AG puede representarse con el siguiente esquema general. Algorithm 6 Algoritmo genético (esquema general) Require: Espacio de búsqueda Ω, aptitud F : Ω→R, tamaño N, probabilidades pc, pm, fracción élite ρ, criterio de paro.
1: Inicializar población P0 = {z(i)}N i=1 ⊂Ωy evaluar F.
2: while no se cumple el criterio de paro do 3:
Seleccionar progenitores (p. ej., torneo) desde P.
4:
Recombinar con prob. pc y mutar con prob. pm; reparar/proyectar si es necesario. 5:
Evaluar la descendencia y formar la nueva población aplicando elitismo (⌊ρN⌋ mejores).
6: end while 7: return mejor individuo observado z⋆.
APÉNDICE G. SINTONÍA ÓPTIMA DE PESOS LQR CON ALGORITMO
GENÉTICO
G.1.
Modelo de optimización Con este andamiaje, el problema se fija en tres piezas: variables de decisión, función objetivo y restricciones.
Variables de decisión.
Se parametrizan los pesos diagonales en escala logarítmica (base 10) para mejorar el acondicionamiento numérico y garantizar positividad: θ = h log10 Q11, . . . , log10 Qnn, log10 R11, . . . , log10 Rmm i⊤∈[ℓb, ub], Q = diag 10θ1, 10θ2, . . . , 10θn
,
R = diag 10θn+1, 10θn+2, . . . , 10θn+m
.
(G.1) Función objetivo.
En el lazo discreto con K = dlqr(Ad, Bd, Q, R) y ganancia de referencia M0 = Cd I −(Ad−BdK) −1Bd, Kr = M †
0,
se minimiza sobre k = 0, . . . , N −1:
m´ın θ J(θ) = Ts X k ∥y[k] −r[k]∥2
2 + λu∥u[k]∥2
2
,
u[k] = −K x[k] + Kr r[k].
Restricciones.
Se aplican tres reglas sencillas: (i) cotas en θ para evitar valores extremos; (ii) positividad de Q, R (implícita al trabajar en escala logarítmica, pues Qii = 10θi y Rjj = 10θn+j); y (iii) estabilidad en lazo cerrado, exigiendo que todos los polos de Ad −BdK queden dentro del círculo unidad. Los individuos que violan (iii) se descartan o se penalizan en J.
θ ∈[ℓb, ub], Q ≻0, R ≻0, ρ(Ad −BdK) < 1.
Acto seguido, se particulariza al caso de estudio LQR discreto con ganancia de referencia Kr, manteniendo un banco de pruebas común para todas las evaluaciones y semillas fijas para reproducibilidad.
G.1.1.
Diagrama de flujo Para clarificar el flujo operativo, la Figura G.1 sintetiza el ciclo evaluación–selección–variación con elitismo y reparación de factibilidad.
Parámetros de entrada (ℓb, ub), N, G, pc, pm, ρ Inicio Inicializar población θ ∼U(ℓb, ub) Evaluar costos J (Estabilidad, K, Kr, simulación) Tomar élite ρN (incumbente) Cruza (prob. pc) uniforme / 1 punto Mutación (prob. pm) gaussiana acotada Evaluar ⇒reparar / proyectar ¿Criterio de paro?
Fin y mejor Q⋆, R⋆ Sí No Figura G.1. Flujo del AG para sintonía de Q, R en LQR. Parámetros y rangos La Tabla G.1 resume la configuración empleada. Los rangos en escala logarítmica para θ ∈[−4, 3] implican Qii, Rjj ∈[10−4, 103], cubriendo sintonías desde suaves hasta agresivas sin deteriorar la condición numérica de la Riccati.
APÉNDICE G. SINTONÍA ÓPTIMA DE PESOS LQR CON ALGORITMO
GENÉTICO
Tabla G.1. Parámetros del algoritmo genético utilizados en los experimentos. Parámetro Valor Comentario Tamaño de población N = 50 Diversidad razonable Generaciones máximas G = 60 Presupuesto de cómputo Élite ρ = 0,2 (Nélite = 10) Preserva incumbente Cruza pc = 0,8 Recombinación uniforme Mutación pm = 0,15 Mutación gaussiana Desviación de mutación σ = 0,25 En escala logarítmica Rangos θ [−4, 3] Qii, Rjj ∈[10−4, 103] Semilla
42
Reproducibilidad Convergencia del costo La evolución del mejor costo por generación se sintetiza en la Tabla G.2: a partir de la generación ≈26 se observa una meseta estable, indicativa de exploración suficiente bajo las cotas impuestas.
Tabla G.2. Convergencia del mejor costo J por generación en LQR-AG+Kr. Gen Mejor J Gen Mejor J Gen Mejor J Gen Mejor J
1
3.506575e+01
8
2.834327e+01
16
2.783841e+01
24
2.776414e+01
2
2.988099e+01
9
2.822356e+01
17
2.781338e+01
26
2.775931e+01
3
2.907251e+01
10
2.816510e+01
18
2.780055e+01
28
2.775776e+01
4
2.882524e+01
11
2.816497e+01
19
2.779671e+01
40
2.774880e+01
5
2.856110e+01
14
2.795018e+01
21
2.777611e+01
43
2.774046e+01
6
2.844072e+01
15
2.788055e+01
23
2.776453e+01
50
2.773273e+01
7
2.836707e+01
13
2.811779e+01
25
2.776363e+01
60
2.773273e+01 Seguimiento de referencia Con K y Kr óptimos fijados, se valida el seguimiento multi-step en las seis salidas. Las Figuras G.2 y G.3 muestran (x, y, z) y (ϕ, θ, ψ), respectivamente, con nudos de cambio marcados para facilitar la lectura de transitorios. La respuesta exhibe transitorios acotados y asentamientos acordes con la sintonía; (ϕ, θ) requieren mayor esfuerzo, sin comprometer estabilidad.
LQR sintonizado por AG (multi-step): (x; y; z) vs referencia
0
5
10
15
20
25
30
Tiempo [s]
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
Posicion [m] x y z Figura G.2. LQR sintonizado por AG (multi–step): (x, y, z) vs referencia. Métricas La Tabla G.3 resume el desempeño del controlador, incluyendo métricas globales y por canal del error y −r.
Tabla G.3. Métricas de desempeño globales y por canal con LQR-AG+Kr. Salida
IAE
ISE
ITAE
RMSE
Eu x
4.1679
3.5017
57.003
0.3163
– y
5.2329
8.6454
104.44
0.4970
– z
0.7915
0.3709
6.439
0.1029
– ψ
1.3984
1.4763
27.533
0.2054
– Global / Total
11.591
13.994
195.42
0.2804
28.498
Los resultados evidencian un desempeño competitivo en términos de seguimiento. El ITAE indica una respuesta adecuada ante cambios en la referencia, mientras que la energía de control Eu refleja el esfuerzo requerido por la sintonía optimizada mediante AG.
APÉNDICE G. SINTONÍA ÓPTIMA DE PESOS LQR CON ALGORITMO
GENÉTICO
LQR sintonizado por AG (multi-step): (?; 3; A) vs referencia
0
5
10
15
20
25
30
Tiempo [s]
-40
-20
0
20
40
60
Angulo [/] ?
3
A Figura G.3. LQR sintonizado por AG (multi–step): (ϕ, θ, ψ) vs referencia.
Apéndice H Control Lineal Cuadrático Gaussiano
(LQG)
El esquema LQG ofrece una vía práctica para operar con sensores ruidosos y estados no medibles directamente: combina una ley de control LQR con un estimador de Kalman estacionario en tiempo discreto. La idea es modular: el LQR fija el compromiso desempeño–energía y el LQE reconstruye estados a partir de salidas afectadas por ruido. Bajo hipótesis gaussianas y condiciones estándar de estabilizabilidad/detectabilidad, el principio de separación garantiza estabilidad de la interconexión. Con este marco, trabajaremos en tiempo discreto con el modelo estocástico xk+1 = Adxk + Bduk + wk, yk = Cdxk + Dduk + vk, donde wk ∼N(0, Qn) y vk ∼N(0, Rn). El desempeño se evalúa mediante J = E " ∞ X k=0 x⊤ k Qcxk + u⊤ k Rcuk #
,
Qc ⪰0, Rc ≻0.
La política adopta uk = −K ˆxk, con K = dlqr(Ad, Bd, Qc, Rc), ˆxk+1 = Adˆxk + Bduk + L yk −Cdˆxk −Dduk
,
y L = dlqe(Ad, I, Cd, Qn, Rn). Si Ad −BdK y Ad −LCd son estables, la interconexión también lo es por separación. Para referencias constantes se emplea una ganancia de referencia estática unitaria, M0 ≜Cd I −(Ad −BdK) −1Bd, Kr = M †
0,
que asegura error nulo en régimen en los canales seleccionados. Con K (realimentación), L (observador) y Kr (ganancia de referencia) queda definida la dinámica de lazo cerrado y, por tanto, la respuesta de salida ante cualquier referencia.
APÉNDICE H. CONTROL LINEAL CUADRÁTICO GAUSSIANO (LQG)
H.1.
Resultados Se consideran dos ensayos complementarios con la misma ganancia de referencia y observador estacionario: primero, un escalón unitario en los canales de posición, manteniendo las referencias angulares en cero; luego, una referencia multi–step por tramos. Así se observa, en orden, la respuesta elemental y la capacidad de seguimiento ante cambios secuenciados.
1) Ensayo a escalón. La Figura H.1 ilustra la respuesta conjunta: (x, y, z) muestran tran-
sitorios moderados y asentamiento consistente con la sintonía; (ϕ, θ, ψ) se mantienen cercanos a cero, por lo que, frente a una consigna unitaria en esos canales, el error final aparece próximo a −1, conforme a la ponderación del costo.
LQG: Seguimiento de Referencia (escalón)
0
2
4
6
8
10
12
14
Tiempo [s]
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
1.2
Posicion [m] / Angulo [rad] x rx y ry z rz ?
r?
3
r3 A rA Figura H.1. LQG: escalón unitario en seis salidas. Transitorios acotados en (x, y, z) y variables angulares próximas a cero por diseño del costo.
Para cuantificar lo observado, la Tabla H.1 resume el desempeño del ensayo a escalón mediante métricas globales y por canal del error y −r.
Tabla H.1. Métricas de desempeño globales y por canal con LQG+Kr (escalón). Salida
IAE
ISE
ITAE
RMSE
Eu x
0.9934
0.7733
0.5779
0.2349
– y
0.9459
0.7369
0.5200
0.2294
– z
3.3047
2.1851
8.0790
0.3949
– ψ 9,5007×10−14 3,4370×10−27 2,1920×10−13 1,5663×10−14 – Global / Total
5.244
3.6953
9.1769
0.2148
0.0664
2) Ensayo multi–step. La Figura H.2 recoge el seguimiento ante cambios por tramos en
(x, y, z), manteniendo ϕ, θ, ψ ancladas en cero. Los nudos verticales marcan los instantes de cambio; la respuesta permanece estable y alineada con el ajuste previo. LQG: Seguimiento de Referencia (multi-step)
0
5
10
15
20
25
30
Tiempo [s]
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
Posicion [m] / Angulo [rad] x rx y ry z rz ?
r?
3
r3 A rA Figura H.2. LQG: seguimiento multi–step. Los ángulos se conservan cercanos a cero y sólo intervienen transitoriamente para sostener la traslación. Las métricas de este experimento se resumen en la Tabla H.2, manteniendo el mismo formato del ensayo a escalón para facilitar la comparación.
APÉNDICE H. CONTROL LINEAL CUADRÁTICO GAUSSIANO (LQG)
Tabla H.2. Métricas de desempeño globales y por canal con LQG+Kr (multi–step). Salida
IAE
ISE
ITAE
RMSE
Eu x
3.9881
3.4734
54.352
0.3150
– y
4.7257
8.1403
93.318
0.4823
– z
5.5578
2.6853
57.610
0.2770
– ψ
1.8964
1.8505
36.714
0.2299
– Global / Total
16.168
16.150
242.00
0.3261
0.1900
Las figuras y tablas muestran estabilidad interna y seguimiento acotado bajo ambos ensayos. En particular, las métricas globales permiten comparar directamente el LQG+Kr con los demás controladores evaluados, mientras que las métricas por canal identifican la contribución específica de x, y, z y ψ al error total.
Bibliografía [1] MarketsandMarkets, “Uav (drone) market.”
https://www.
marketsandmarkets.com, 2024.
[2] Drone Industry Insights, “Global drone market report 2023-2030,” MarketsandMarkets, 2023.
[3] NASA, “NASA’s Drone Market Report.” https://www.nasa.gov, 2023. [4] Naciones Unidas, “Objetivos de desarrollo sostenible (ods).”
https://www.un.org/sustainabledevelopment/es/ objetivos-de-desarrollo-sostenible/, 2024. [5] Federación Internacional de Robótica (IFR), “Informe sobre la adopción de drones en américa latina.” https://ifr.org, 2022. [6] Cámara Colombiana de Comercio Electrónico (CCCE), “Informe de crecimiento del uso de drones en colombia.” https://www.ccce.org.co, 2021. [7] Departamento Nacional de Planeación, “Plan nacional de desarrollo 2022-2026.” https://www.dnp.gov.co, 2023. ISSN: 2022-2026. [8] Statista Market Insights, “Drones: Market data analysis.” https://www.
statista.com, 2024.
[9] NASA,
“About pathfinding for airspace with autonomous vehicles.”
https://www.nasa.gov/directorates/armd/aosp/atm-x/paav/ about-paav/, 2024.
[10] S. Nahavandi, R. Alizadehsani, D. Nahavandi, S. Mohamed, N. Mohajer, M. Rokonuzzaman, and I. Hossain, “A comprehensive review on autonomous navigation,” ACM Comput. Surv., vol. 57, May 2025.
BIBLIOGRAFÍA
[11] J. Van Brummelen, M. O’Brien, D. Gruyer, and H. Najjaran, “Autonomous vehicle perception: The technology of today and tomorrow,” Transportation Research Part C: Emerging Technologies, vol. 89, pp. 384–406, 2018.
[12] S. Zhang, Y. Li, and Q. Dong, “Autonomous navigation of uav in multi-obstacle environments based on a deep reinforcement learning approach,” Applied Soft Computing, vol. 115, p. 108194, 2022.
[13] B. Lindqvist, S. Karlsson, A. Koval, I. Tevetzidis, J. Haluška, C. Kanellakis, A.-A. Agha-Mohammadi, and G. Nikolakopoulos, “Multimodality robotic systems: Integrated combined legged-aerial mobility for subterranean search-and-rescue,” Robotics and Autonomous Systems, vol. 154, p. 104134, 2022.
[14] C. Cheng, Q. Sha, B. He, and G. Li, “Path planning and obstacle avoidance for auv: A review,” Ocean Engineering, vol. 235, 9 2021.
[15] V. Yordanov, L. Barazzetti, M. A. Brovelli, J. Sun, G. Yuan, L. Song, and H. Zhang, “Unmanned aerial vehicles (uavs) in landslide investigation and monitoring: A review,” Drones 2024, Vol. 8, Page 30, vol. 8, p. 30, 1 2024. [16] A. A. Laghari, A. K. Jumani, R. A. Laghari, and H. Nawaz, “Unmanned aerial vehicles: A review,” Cognitive Robotics, vol. 3, pp. 8–22, 2023. [17] D. Falanga, K. Kleber, S. Mintchev, D. Floreano, and D. Scaramuzza, “The foldable drone: A morphing quadrotor that can squeeze and fly,” IEEE Robotics and Automation Letters, vol. 4, no. 2, pp. 209–216, 2019.
[18] H.-Y. Lin and X.-Z. Peng, “Autonomous quadrotor navigation with vision based obstacle avoidance and path planning,” IEEE Access, vol. 9, pp. 102450–102459,
2021.
[19] T. GUO, N. JIANG, B. LI, X. ZHU, Y. WANG, and W. DU, “Uav navigation in high dynamic environments: A deep reinforcement learning approach,” Chinese Journal of Aeronautics, vol. 34, no. 2, pp. 479–489, 2021.
[20] A. Romero, R. Penicka, and D. Scaramuzza, “Time-optimal online replanning for agile quadrotor flight,” IEEE Robotics and Automation Letters, vol. 7, pp. 7730–
7737, 7 2022.
[21] A. Marashian and A. Razminia, “Mobile robot’s path-planning and path-tracking in static and dynamic environments: Dynamic programming approach,” Robotics and Autonomous Systems, vol. 172, p. 104592, 2024.
[22] D. Falanga, K. Kleber, and D. Scaramuzza, “Dynamic obstacle avoidance for quadrotors with event cameras,” Science Robotics, vol. 5, p. eaaz9712, 2020. [23] C. Yin, Z. Xiao, X. Cao, X. Xi, P. Yang, and D. Wu, “Offline and online search: Uav multiobjective path planning under dynamic urban environment,” IEEE Internet of Things Journal, vol. 5, pp. 546–558, 4 2018.
[24] J. Li, X. Xiong, Y. Yan, and Y. Yang, “A survey of indoor uav obstacle avoidance research,” IEEE Access, vol. 11, pp. 51861–51891, 2023.
[25] J. N. Yasin, S. A. S. Mohamed, M.-H. Haghbayan, J. Heikkonen, H. Tenhunen, and J. Plosila, “Unmanned aerial vehicles (uavs): Collision avoidance systems and approaches,” IEEE Access, vol. 8, pp. 105139–105155, 2020.
[26] B. Zhou, J. Pan, F. Gao, and S. Shen, “Raptor: Robust and perception-aware trajectory replanning for quadrotor fast flight,” IEEE Transactions on Robotics, vol. 37, pp. 1992–2009, 2021.
[27] M. Y. Arafat, M. M. Alam, and S. Moh, “Vision-based navigation techniques for unmanned aerial vehicles: Review and challenges,” Drones, vol. 7, 2 2023. [28] F. Gao, L. Wang, B. Zhou, X. Zhou, J. Pan, and S. Shen, “Teach-repeat-replan: A complete and robust system for aggressive flight in complex environments,” IEEE Transactions on Robotics, vol. 36, pp. 1526–1545, 10 2020.
[29] J. Wu, Y. Ye, and J. Du, “Multi-objective reinforcement learning for autonomous drone navigation in urban areas with wind zones,” Automation in Construction, vol. 158, 2 2024.
[30] P. S. Fakhri, O. Asghari, S. Sarspy, M. B. Marand, P. Moshaver, and M. Trik, “A fuzzy decision-making system for video tracking with multiple objects in nonstationary conditions,” Heliyon, vol. 9, p. e22156, November 2023. [31] M. Z. Butt, N. Nasir, R. Bt, and A. Rashid, “A review of perception sensors, techniques, and hardware architectures for autonomous low-altitude uavs in noncooperative local obstacle avoidance,” Robotics and Autonomous Systems, vol. 173,
p. 104629, 2024.
BIBLIOGRAFÍA
[32] G. Godinez-Garrido and J. Gonzalez-Islas, “Estimation of damaged regions by the bark beetle in a mexican forest using uav images and deep learning,” Sustainability, vol. 16, no. 23, p. 10731, 2024.
[33] M. Ahmed, M. Adnan, M. Ahmed, and D. Janssens, “From stationary to nonstationary uavs: Deep-learning-based method for vehicle speed estimation,” Algorithms, vol. 17, no. 12, p. 558, 2024.
[34] S. Han and D. Han, “Enhancing direct georeferencing using real-time kinematic uavs and structure from motion-based photogrammetry for large-scale infrastructure,” Drones, vol. 8, no. 12, p. 736, 2024.
[35] A. Jamali, B. Lu, E. M. Gerbrandt, C. Teasdale, R. R. Burlakoti, S. Sabaratnam, J. McIntyre, L. Yang, M. Schmidt, D. McCaffrey, and P. Ghamisi, “High-resolution uav-based blueberry scorch virus mapping utilizing a deep vision transformer algorithm,” Computers and Electronics in Agriculture, vol. 229, p. 109726, 2025. [36] Y. Wang, W. Zhao, R. Zhang, N. Li, D. Li, and J. Lv, “Multi-object tracking in uavs with feature fusion distribution and occlusion awareness,” Signal, Image and Video Processing, vol. 19, pp. 453–469, 2025.
[37] H. Chen, K. Wu, H. Lin, H. Zhou, Z. Zhou, and Y. Mai, “A real-time vision guidance method for autonomous longan picking by the uav,” Computers and Electronics in Agriculture, vol. 211, pp. 128–145, 2025.
[38] Y. Chang, Y. Cheng, U. Manzoor, and J. Murray, “A review of uav autonomous navigation in gps-denied environments,” Robotics and Autonomous Systems, vol. 170,
p. 104533, 2023.
[39] K. Pereida and A. Schoellig, “Adaptive model predictive control for high-accuracy trajectory tracking in changing conditions,” IEEE Transactions on Robotics, vol. 40, no. 1, pp. 12–25, 2024.
[40] W. Zheng and B. Zhu, “Stochastic time-varying model predictive control for trajectory tracking of a wheeled mobile robot,” Frontiers in Energy Research, vol. 11,
p. 215, 2023.
[41] R. Kumar, P. Verma, and L. Gupta, “Energy-efficient trajectory planning for uav networks using deep reinforcement learning,” IEEE Access, vol. 12, pp. 20345–20358,
2024.
[42] V. Sankaranarayanan, G. Damigos, A. S. Seisa, S. Satpute, T. Lindgren, and G. Nikolakopoulos, “Paced-5g: Predictive autonomous control using edge for drones over 5g,” IEEE Internet of Things Journal, vol. 10, no. 2, pp. 89–101, 2023. [43] Y. Wang, J. O’Keeffe, Q. Qian, and D. Boyle, “Interpretable stochastic model predictive control using distributional reinforced estimation for quadrotor tracking systems,” IEEE Transactions on Cybernetics, vol. 53, no. 4, pp. 562–580, 2022. [44] T. Miller and S. Johnson, “Anticipatory path planning for multi-uav systems using predictive analytics,” Journal of Intelligent Robotic Systems, vol. 95, pp. 305–320,
2024.
[45] M. Kazim, H. Sim, G. Shin, H. Hwang, and K. Kim, “Aggressive trajectory tracking for nano quadrotors using embedded nonlinear model predictive control,” IEEE Robotics and Automation Letters, vol. 8, no. 1, pp. 101–115, 2023. [46] C. Li and F. Zhang, “Model-free adaptive control for uav swarms in dynamic environments,” IEEE Transactions on Control Systems Technology, vol. 33, no. 4, pp. 4021–4036, 2025.
[47] D. Hanover, A. Loquercio, L. Bauersfeld, A. Romero, R. Penicka, Y. Song, G. Cioffi, E. Kaufmann, and D. Scaramuzza, “Autonomous drone racing: A survey,” Trans. Rob., vol. 40, p. 3044–3067, Jan. 2024.
[48] A. Soler, J. Betancurt, and D. Amortegui, “Diseño e implementación de un controlador multivariable en un dispositivo UAV tipo cuadricóptero para fumigación aérea,” trabajo de grado, Universidad Tecnológica de Pereira, Pereira, Colombia, 2016. [49] J. A. B. Becerra, “Control robusto acoplado de sistemas multivariables aplicado a vehículos aéreos no tripulados,” trabajo de grado, Universidad Tecnológica de Pereira, Pereira, Colombia, 2019.
[50] A. A. M. Lopez, Control multivariable no lineal de estructura variable aplicado a un vehículo aéreo no tripulado. PhD thesis, Universidad Tecnológica de Pereira, Pereira, Colombia, 2017.
[51] M. L. Rivera, “Implementación de técnicas de control inteligente en un helicóptero no tripulado de dos grados de libertad,” trabajo de grado, Universidad Tecnológica de Pereira, Pereira, Colombia, 2020.
BIBLIOGRAFÍA
[52] Y. Liao, Y. Wu, S. Zhao, and D. Zhang, “Unmanned aerial vehicle obstacle avoidance based custom elliptic domain,” Drones, vol. 8, no. 8, 2024.
[53] S. Hutchinson, G. D. Hager, and P. I. Corke, “A tutorial on visual servo control,” IEEE Transactions on Robotics and Automation, vol. 12, no. 5, pp. 651–670, 1996. [54] F. Chaumette and S. Hutchinson, “Visual servo control. i. basic approaches,” IEEE Robotics & Automation Magazine, vol. 13, no. 4, pp. 82–90, 2006. [55] F. Chaumette and S. Hutchinson, “Visual servo control. ii. advanced approaches [tutorial],” IEEE Robotics & Automation Magazine, vol. 14, no. 1, pp. 109–118, 2007. [56] A. Keipour, G. A. S. Pereira, R. Bonatti, R. Garg, P. Rastogi, G. Dubey, and S. Scherer, “Visual servoing approach to autonomous UAV landing on a moving vehicle,” Sensors, vol. 22, no. 17, p. 6549, 2022.
[57] G. Cho, S.-H. Choi, J. Bae, and H. Oh, “Autonomous ship deck landing of a quadrotor UAV using feed-forward image-based visual servoing,” Aerospace Science and Technology, vol. 130, p. 107869, 2022.
[58] J. Wu, Z. Jin, A. Liu, L. Yu, and F. Yang, “A survey of learning-based control of robotic visual servoing systems,” Journal of the Franklin Institute, vol. 359, no. 1, pp. 556–577, 2022.
[59] F. Prochazka, S. Krüger, G. Stomberg, and M. Bauer, “Development of a hardwarein-the-loop demonstrator for the validation of fault-tolerant control methods for a hybrid uav,” CEAS Aeronautical Journal, vol. 12, no. 2, pp. 123–135, 2021. [60] S. Park, H. Kim, and H. Kim, “Real-time validation of formation control for fixedwing uavs using hardware-in-the-loop simulation,” IEEE Access, vol. 7, pp. 34562–
34575, 2019.
[61] Y. Zhang, X. Li, J. Wang, and L. Zhang, “Hardware-in-the-loop simulation platform for unmanned aerial vehicle swarm system: Architecture and application,” IEEE Transactions on Industrial Informatics, vol. 16, no. 3, pp. 567–580, 2020. [62] M. Mammarella and E. Capello, “Tube-based robust mpc processor-in-the-loop validation for fixed-wing uavs,” Journal of Intelligent Robotic Systems, vol. 98, pp. 123–
140, 2020.
[63] J. Liu, Y. Wang, and H. Zhang, “Hardware-in-the-loop based 6dof test platform for multi-rotor uav,” IEEE Transactions on Aerospace and Electronic Systems, vol. 55, no. 4, pp. 3456–3470, 2019.
[64] A. Smith, B. Johnson, and C. Lee, “Comprehensive hardware-in-the-loop simulation architecture for quadrotor helicopters,” IEEE Transactions on Control Systems Technology, vol. 30, no. 5, pp. 3500–3515, 2022.
[65] Z. Xu, X. Zhan, B. Chen, Y. Xiu, C. Yang, and K. Shimada, “A real-time dynamic obstacle tracking and mapping system for uav navigation and collision avoidance with an rgb-d camera,” in 2023 IEEE International Conference on Robotics and Automation (ICRA), pp. 10645–10651, IEEE, 2023.
Cita: Quintana-Fuentes, Jose Daniel (2026), Sistema integrado de percepción visual y control adaptativo para navegación autónoma de Quadcopters en entornos dinámicos, Universidad Tecnológica de Pereira, p. N. https://hdl.handle.net/11059/16814