Estudio de integración de sensores en UAVs

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo de Fin de Grado Ingeniería Aeroespacial Estudio de integración de sensores en UAVs Autor: José Luis Carretero Rodríguez Tutora: Juana María Martínez Heredia Departamento de Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 Trabajo de Fin de Grado Ingeniería Aeroespacial Estudio de integración de sensores en UAVs Autor: José Luis Carretero Rodríguez Tutora: Juana María Martínez Heredia Departamento de Ingeniería Electrónica Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 Trabajo de Fin de Grado: Estudio de integración de sensores en UAVs Autor: José Luis Carretero Rodríguez Tutora: Juana María Martínez Heredia El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Índice general 1. Introducción 4 2. Modelado en el Espacio de Estados 7 2.1. Representación de sistemas en el espacio estados . . . . . . . . . . . . . . . . . . . . . 8 2.2. Obtención de la representación en espacio de estados de sistemas discretos . . . . . . . 9 2.2.1. Método de programación directa . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2. Método de programación anidada . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3. Relación entre la representación en espacio de estados y la función de transferencia . . 13 2.4. No unicidad de la representación en espacio de estados de un sistema . . . . . . . . . . 15 2.5. Resolución de las ecuaciones del espacio de estados . . . . . . . . . . . . . . . . . . . . 15 2.5.1. Procedimiento recursivo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5.2. Matriz de transición de estados . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5.3. Método basado en la transformada Z . . . . . . . . . . . . . . . . . . . . . . . . 17 2.6. Linealización de las ecuaciones de estado . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.6.1. Interpretación Gráfica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.6.2. Aproximación lineal del modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.7. Discretización de las ecuaciones de estado continuas . . . . . . . . . . . . . . . . . . . . 20 2.8. Controlabilidad y Observabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.8.1. Controlabilidad del estado completo . . . . . . . . . . . . . . . . . . . . . . . . 25 2.8.2. Controlabilidad de la salida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.8.3. Observabilidad .................................... 27 2.8.4. Principio de Dualidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.9. Transformación de un sistema en formas canónicas . . . . . . . . . . . . . . . . . . . . 29 2.9.1. Forma Canónica Controlable . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.9.2. Forma Canónica Observable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.10. Descripción de un sistema en parte controlable/observable y no controlable/no observable 31 2.10.1. Parte controlable/no controlable . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.10.2. Parte observable/no observable . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3. Modelado de Sensores en los Vehículos Aéreos no Tripulados 34 3.1. Caracterización de los Sensores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 ii 0 ÍNDICE GENERAL 3.1.1. Descriptores estáticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1.2. Descriptores dinámicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2. Uso correcto de los sensores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.1. Calibración de sensores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2.2. Modelado de simulaciones SIL-HIL - Software-in-the-loop yHardware-in-the-loop 39 3.3. Sensores para medir distancias y proximidad . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.1. Sensores capacitivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.2. Sensoresinductivos.................................. 42 3.3.3. Basados en efecto Hall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.4. Basados en ultrasonidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3.5. Sensores de Espectro Infrarrojo . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3.6. Visiónartificial.................................... 45 3.4. Sensoresdeluz ........................................ 47 3.5. SensoresdeVelocidad .................................... 48 3.5.1. TubodePitot..................................... 50 3.6. Codificadores (Encoders)................................... 51 3.7. Altímetros........................................... 52 3.8. Giroscopios, Acelerómetros y Magnetómetros: Los pilares de una IMU (del inglés, Inertial Measurement Unit - Unidad Inercial de Medida) . . . . . . . . . . . . . . . . . . . 52 3.9. GPS - Global Positioning System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.9.1. Segmentos....................................... 57 3.9.2. Disponibilidad, Integridad y Continuidad . . . . . . . . . . . . . . . . . . . . . 58 4. Filtro de Kalman 60 4.1. Introducción.......................................... 60 4.2. ¿QuéeselfiltrodeKalman?................................. 60 4.3. El modelo multidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5. Caso práctico - Simulación en Matlab del Filtro de Kalman 68 5.1. Procesamientodedatos ................................... 68 5.1.1. Medida de la trayectoria - GPS . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.1.2. Medida de la trayectoria - Google Maps . . . . . . . . . . . . . . . . . . . . . . 70 5.2. Parámetros del filtro de Kalman . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.3. Resultados obtenidos de la simulación . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 6. Conclusiones y Líneas Futuras 76 A. Varianza & Covarianza 78 B. Código Matlab de caso práctico 79 iii 0 ÍNDICE GENERAL Bibliografía 86 1 2 Modelado en el Espacio de Estados 2.1. Representación de sistemas en el espacio estados Con la representación en espacio de estados se obtiene la capacidad de conocer y controlar en cierta medida la dinámica interna de un sistema y su respuesta. Las variables contenidas en el vector de estado deben ser capaces en su conjunto de determinar las condiciones de la dinámica del sistema. Cabe destacar en este punto que es posible que existan varias representaciones en variables de estado para un mismo sistema. Para una mejor comprensión de los conceptos que aquí se tratarán, se procede a definir brevemente los siguientes términos [4], [13]: Sistema, se entenderá como una relación entre entradas y salidas. Sistema Determinista, si a cada entrada le corresponde una y solo una salida Sistema Monovariable. Es aquel que solo tiene una entrada y una salida. Si el sistema tiene más de una entrada o más de una salida se le considerará multivariable. Sistema Causal o No Anticipatorio. Es aquel que su salida para cierto tiempo t1, no depende de entradas aplicadas después de t1. Obsérvese que la definición implica que un sistema no causal es capaz de predecir entradas futuras, por lo tanto la causalidad es una propiedad intrínseca de cualquier sistema físico. Sistema Dinámico. Es aquel cuya salida presente depende de entradas pasadas y presentes. Si el valor de la salida en t1depende solamente de la entrada aplicada en t1, el sistema se conoce como estático o sin memoria. La salida de un sistema estático permanece constante si la entrada no cambia. En un sistema dinámico la salida cambia con el tiempo aunque no se cambie la entrada, a menos que el sistema ya se encuentre en estado estable. Sistema Invariante en el Tiempo. Es aquel que tiene parámetros fijos o estacionarios con respecto al tiempo, es decir, sus características no cambian al pasar el tiempo o dicho de otra forma, sus propiedades son invariantes con traslaciones de tiempo. Sistema No Lineal. Un sistema es no lineal si no se aplica el principio de superposición. Por ello, para un sistema que sea no lineal, la respuesta a la suma de dos entradas no puede calcularse tratando por separado a cada entrada y sumando los resultados obtenidos. La dinámica de un sistema se puede describir en función del valor del vector de estados y de la señal de entrada (asumiendo que el sistema es no autónomo) mediante las siguientes expresiones: x(k+ 1) = f(x(k), u(k), k) y(k) = g(x(k), u(k), k)(2.2) donde la notación ξkindica el valor tomado por ξen el instante de tiempo tkyfygpueden ser funciones de cualquier tipo. Otra posible representación de este tipo de sistemas es la siguiente: x(k+ 1) = Gx(k) + Hu(k) y(k) = Cx(k) + Du(k)(2.3) la cual es utilizada para trabajar con sistemas lineales e invariantes en el tiempo (LTI), y donde x(k) = Vector de estado (dimensión n) y(k) = Vector de salida (dimensión m) 8 2 Modelado en el Espacio de Estados u(k) = Vector de entrada (dimensión r) G(k) = Matriz de estado (dimensión n×n) H(k) = Matriz de entrada (dimensión n×r) C(k) = Matriz de salida (dimensión m×n) D(k) = Matriz de transmisión directa (dimensión m×r) Las ecuaciones 2.3 también se pueden representar mediante el diagrama de bloques de la Figura 2.1 Figura 2.1: Diagrama de bloques de la representación en espacio estados de un sistema LTI [13] 2.2. Obtención de la representación en espacio de estados de sistemas discretos Para comenzar esta sección, se explicará previamente la diferencia que existe entre un sistema continuo y un sistema discreto. Sistema en tiempo Continuo. Es aquel en el que las señales con las que trabaja el sistema son funciones de la variable continua tiempo t. Un ejemplo de señal continua [x(t)] sería la variación de temperatura que se produce en los componentes de un sistema, la intensidad luminosa que se recibe del sol, señales sinusoidales, etc. Sistema en tiempo Discreto. Es aquel sistema en el cual una o más de las variables involucradas pueden cambiar solamente en valores concretos de tiempo. A dichos instantes se les denotará como kokT, e indican los tiempos en los que se lleva a cabo alguna medición de tipo físico o el instante en el que se extraen datos de algún parámetro, variable u otro sistema. Es decir, una señal discreta [x(kT)] solo toma valores para ciertos puntos, de forma que entre dos instantes en los que se toman valores se mantendrá el último valor tomado. Las señales digitales son el mejor ejemplo de señales discretas. Como ejemplos, las señales de GPS son discretas porque se envían cada cierto periodo de tiempo, los pulsos cardíacos y los semáforos, entre otras. En la sección 2.7 se explicará el proceso mediante el cual una señal continua podrá ser discretizada para su uso dentro del modelado en espacio estados. En el caso bajo estudio, de la simulación de la integración de dos señales de GPS y de IMU de un UAV, se utilizará un sistema en tiempo discreto. Esto será debido a que las señales tanto por parte del GPS como de la IMU se reciben en instantes de tiempo concretos, discretos. Se puede decir que el sistema discreto bajo estudio es descrito por la siguiente expresión: y(k) + a1y(k−1) + a2y(k−2) + ... +any(k−n) = b0u(k) + b1u(k−1) + ... +bnu(k−n)(2.4) 9 2 Modelado en el Espacio de Estados donde u(k) es la entrada (con sus correspondientes coeficientes bj) e y(k) es la salida (con sus correspondientes coeficientes aj) del sistema en el instante de muestreo k. Es posible que alguno de los coeficientes ai(i=1,2,...,n) y bj(j=0,1,2,...,n) pueden ser cero. En otras palabras, esta expresión significa que la salida del sistema actual y(k)(lo que por ejemplo podría ser la posición actual del vehículo), depende tanto de las posiciones anteriores y(k−n)y de las variables de entrada u(k)(las cuales podrían ser una fuerza que desplazara el vehículo como el viento, la gravedad, una variación de velocidad desde los mandos de control, etc.). A partir de su función de transferencia será posible obtener la representación de espacio de estados del sistema. Puesto que la función de transferencia proporciona la relación que existe entre las entradas y salidas del sistema, se obtiene: G(z) = Y(z) U(z)=b0+b1z−1+b2z−2+... +bnz−n 1 + a1z−1+a2z−2+... +anz−n(2.5) Existen dos métodos para obtener la representación de espacio de estados a partir de 2.5: el Método de programación directa y el Método de programación anidada, los cuales se explicarán a continuación. [13] [16] [25]. 2.2.1. Método de programación directa Se reescribe la función de transferencia como: G(z) = b0+(b1−a1b0)z−1+ (b2−a2b0)z−2+... + (bn−anb0)z−n 1 + a1z−1+a2z−2+... +anz−n(2.6) y también teniendo en cuenta que G(z) = Y(z) U(z), se obtiene: Y(z) = b0U(z) + (b1−a1b0)z−1+ (b2−a2b0)z−2+... + (bn−anb0)z−n 1 + a1z−1+a2z−2+... +anz−nU(z)(2.7) que puede reescribirse como: Y(z) = b0U(z) + ˜ Y(z)U(z)⇔˜ Y(z) = (b1−a1b0)z−1+ (b2−a2b0)z−2+... + (bn−anb0)z−n 1 + a1z−1+a2z−2+... +anz−n(2.8) Conociendo la expresión de ˜ Y(z)se puede conseguir una variable auxiliar intermedia entre las entradas y salidas, Q(z), a partir de la cual se podrán definir las relaciones entre el estado en un instante y en el instante anterior (como se verá en la formulación del modelo de espacio estados de la expresión 2.13), y que además cumple lo siguiente Q(z) = ˜ Y(z) (b1−a1b0)z−1+ (b2−a2b0)z−2+... + (bn−anb0)z−n=U(z) 1 + a1z−1+a2z−2+... +anz−n (2.9) De ahí se obtiene que: Q(z) = −a1z−1Q(z)−a2z−2Q(z)−... −anz−nQ(z) + U(z)(2.10) ˜ Y(z)=(b1−a1b0)z−1Q(z)+(b2−a2b0)z−2Q(z) + ... + (bn−anb0)z−nQ(z)(2.11) Las variables de estado del problema serían: 10 2 Modelado en el Espacio de Estados X1(z) = z−nQ(z) X2(z) = z−(n−1)Q(z) ... Xn(z) = z−1Q(z) Utilizando las propiedades de la transformada de Z: zX1(z) = X2(z)⇔x1(k+ 1) = x2(k) zX2(z) = X3(z)⇔x2(k+ 1) = x3(k) ... ⇔... zX(n−1)(z) = Xn(z)⇔x(n−1)(k+ 1) = xn(k) Teniendo en cuenta las expresiones anteriores, es posible obtener Q(z) = zXn(z), y por tanto se puede reescribir la expresión de 2.10 como: zXn(z) = −a1Xn(z)−a2X(n−1)(z)−... −anX1(z) + U(z)⇔ xn(k+ 1) = −anx1(k)−an−1x2(k)−... −a1xn(k) + u(k)(2.12) Con esta información ya es posible obtener aquello que se estaba buscando: la expresión de la ecuación de estado. De forma matricial, quedaría de la siguiente forma:        x1(k+ 1) x2(k+ 1) . . . x(n−1)(k+ 1) xn(k+ 1)        =        0 1 0 ... 0 0 0 1 ... 0 . . .. . .. . ..... . . 0 0 0 ... 1 −an−an−1−an−2... −a1               x1(k) x2(k) . . . x(n−1)(k) xn(k)        +        0 0 . . . 0 1        uk(2.13) En la ecuación 2.13 se puede observar cómo el siguiente estado (x1(k+ 1)) depende del anterior (x2(k)). Esta ecuación es una expresión bastante clara donde se pueden observar los estados internos del sistema, de una forma bastante visual, fenómeno que no se producía cuando se tenía cuando se observaba directamente la función de transferencia. La ecuación 2.13 es la nombrada anteriormente como “Ecuación de estado”. Si del sistema bajo estudio se quieren analizar los datos de la posición y velocidad (estados), por ejemplo, pues esta expresión sería la necesaria para calcular dichos datos, los cuales se basan en los estados anteriores y en las entradas. Para calcular ahora la denominada “Ecuación de salida” procedente de la función de transferencia anterior, simplemente bastaría con reescribir la ecuación 2.11: ˜ Y(z)=(b1−a1b0)Xn(z)+(b2−a2b0)X(n−1)(z) + ... + (bn−anb0)X1(z) ⇓antitransformando y(k)=(bn−anb0)x1(k)+(bn−1−an−1b0)x2(k) + ... + (b1−a1b0)xn(k) + b0u(k) (2.14) Se llega a una expresión de la ecuación de la salida de la siguiente forma: 11 2 Modelado en el Espacio de Estados y(k) = bn−anb0bn−1−an−1b0· · · b1−a1b0        x1(k) x2(k) . . . x(n−1)(k) xn(k)        +b0u(k)(2.15) Como se ha dicho anteriormente al principio de esta sección, existen varias representaciones de espacio de estados para un mismo sistema. Las ecuaciones de estado2.13 y de salida 2.15 son un ejemplo de una de esas representaciones. Ambas representan el espacio estados del sistema cuya función de transferencia es la expresión 2.5, que se denomina forma canónica controlable (FCC). Más adelante se explicarán con más detalle este tipo de representación. 2.2.2. Método de programación anidada Se parte en este caso de la misma función de transferencia 2.5, para encontrar otro tipo de representación del espacio de estados. Y(z)−b0U(z) + z−1(a1Y(z)−b1U(z)) + ... +z−n(anY(z)−bnU(z)) = 0 ⇔ Y(z) = b0U(z) + z−1(b1U(z)−a1Y(z) + z−1(b2U(z)−a2Y(z) + z−1(b3U(z)−a3Y(z) + ...))) ⇔ Y(z) = b0U(z) + Xn(z) =⇒antitransformando =⇒y(k) = xn(k) + b0u(k) (2.16) A partir de esta expresión se pueden definir las siguientes variables de estado: Xn(z) = z−1(b1U(z)−a1Y(z) + X(n−1)(z)) X(n−1)(z) = z−1(b2U(z)−a2Y(z) + X(n−2)(z)) ... X2(z) = z−1(bn−1U(z)−an−1Y(z) + X1(z)) X1(z) = z−1(bnU(z)−anY(z)) (2.17) Si se sustituye la última expresión de 2.16 en las variables de estado, y a su vez se multiplica por z en ambos miembros de la igualdad: zXn(z) = X(n−1)(z)−a1Xn(z)+(b1−a1b0)U(z) zX(n−1)(z) = X(n−2)(z)−a2Xn(z)+(b2−a2b0)U(z) ... zX2(z) = X1(z)−an−1Xn(z)+(bn−1−an−1b0)U(z) zX1(z) = −anXn(z)+(bn−anb0)U(z) (2.18) Antitransformando la ecuación 2.18: x1(k+ 1) = −anxn(k)+(bn−anb0)u(k) x2(k+ 1) = x1(k)−an−1xn(k)+(bn−1−an−1b0)u(k) ... x(n−1)(k+ 1) = x(n−2)(k)−a2xn(k)+(b2−a2b0)u(k) xn(k+ 1) = x(n−1)(k)−a1xn(k)+(b1−a1b0)u(k) (2.19) 12 2 Modelado en el Espacio de Estados Teniendo en cuenta la antitransformada que se obtuvo en la ecuación 2.16 y 2.19, se pueden obtener la representación de espacio de estados de el sistema bajo estudio.        x1(k+ 1) x2(k+ 1) . . . x(n−1)(k+ 1) xn(k+ 1)        =        0 0 · · · 0 0 −an 1 0 · · · 0 0 −an−1 . . .. . ..... . .. . .. . . 0 0 · · · 1 0 −a2 0 0 · · · 0 1 −a1               x1(k) x2(k) . . . x(n−1)(k) xn(k)        +        bn−anb0 bn−1−an−1b0 . . . b2−a2b0 b1−a1b0        u(k) y(k) = 0 0 · · · 0 1         x1(k) x2(k) . . . x(n−1)(k) xn(k)        +b0u(k) (2.20) A este tipo de representación, donde en la expresión 2.20 se aprecian la ecuación de estado y de salida, se le denomina forma canónica observable (FCO). Se estudiará más adelante por qué motivos se utilizará esta forma de representación o la anterior. Cuando se realiza un análisis y diseño de un modelo concreto en el dominio de estado, típicamente se transforman las ecuaciones de las que se dispone en alguna de estas formas particulares debido a las ventajas que ofrecen. Por ejemplo, la FCC posee propiedades interesantes que la hacen conveniente para pruebas de controlabilidad, mientras que la FCO se utiliza más para pruebas de observabilidad. Ambos conceptos (controlabilidad y observabilidad) serán introducidos en la sección 2.8. 2.3. Relación entre la representación en espacio de estados y la función de transferencia El modelado y control de sistemas basado en la transformada de Laplace ofrece un enfoque sencillo y de fácil aplicación. Permite analizar sistemas utilizando una serie de reglas algebraicas en lugar de trabajar con ecuaciones diferenciales. Pero no en todas sus formas tiene la misma elegancia. Las funciones de transferencia cuentan con una serie de limitaciones a la hora de describir un sistema [4]: No proporciona información sobre la estructura física del sistema Sólo es válida para sistemas lineales con una entrada y una salida e invariantes en el tiempo. Este es uno de los motivos por los cuales se descarta trabajar con funciones de transferencia y se opta por el espacio estados. En el sistema bajo estudio, en el cual se integran señales de sistemas GPS e IMU, las señales recibidas son discretas (con las que trabaja el modelado en espacio estados) y no continuas (con las que trabaja la transformada de Laplace). Como se verá en el apartado de Sensores, existen múltiples entradas y múltiples salidas (se obtienen señales de altímetros, velocímetros, cámaras, codificadores, giróscopos, acelerómetros, magnetómetros, etc.) que se querrán integrar unas con otras para obtener el mejor resultado posible, lo cual no puede llevarse a cabo con una función de transferencia (trabaja con una única entrada y única salida). No proporciona información de lo que pasa dentro del sistema. Si se quiere obtener un dato intermedio de alguno de los sensores o sistemas que se encuentran operativos en el proceso, es más complejo obtenerlo en el caso de trabajar con este tipo de descripciones. Se necesita que las condiciones iniciales del sistemas sean nulas. De hecho, la función de transferencia de un sistema lineal invariante en el tiempo se define como la transformada de Laplace de 13 2 Modelado en el Espacio de Estados la respuesta al impulso, con todas las condiciones iniciales iguales a cero. En el caso del vehículo bajo estudio, por ejemplificar, se quieren tomar datos de una trayectoria rectilínea con velocidad constante de 8m/s. Es decir, comenzando y terminando con esa velocidad. En este caso la velocidad inicial no es nula, y por tanto, no se podría modelar este proceso mediante la función de transferencia. Otro motivo por el cual se utilizará el modelado en espacio estados, porque tal y como se ha demostrado, no limita la maniobrabilidad de actuaciones que se pueden llevar a cabo No solamente el caso bajo estudio se encuentra fuera del rango de actuación de la función de transferencia, sino que la mayoría de sistemas dinámicos no cumplen con dichos requisitos. Los sistemas reales, por lo general, presentan no linealidades, cuentan con más de una entrada y salida, sus parámetros cambian con el tiempo y sus condiciones iniciales no siempre son cero. Con un sencillo análisis de la representación de espacio de estados se puede obtener la función de transferencia. ˙x=Ax +Bu y=Cx +Du =⇒Laplace =⇒sX(s) = AX(s) + BU(s) Y(s) = CX(s) + DU(s)(2.21) De la ecuación de estado se puede obtener: (sI −A)X(s) = BU(s) =⇒X(s) = (sI −A)−1BU(s)(2.22) Y sustituyéndolo en la ecuación de salida: Y(s) = [C(sI −A)−1B+D]U(s) =⇒G(s) = Y(s) U(s)= [C(sI −A)−1B+D](2.23) Pero lo interesante en este tema no es obtener la función de transferencia a partir de la representación de espacio estados, sino al contrario. Como se ha repetido anteriormente, el proceso de convertir de función de transferencia a espacio de estados no es único. Se puede decir que todas las transformaciones son ”equivalentes”, puesto que las propiedades del sistema no cambian. Lo atractivo de este proceso es que algunas representaciones de espacio de estados pueden tener mas ventajas que otras dependiendo del caso para una tarea particular. Algunas posibles representaciones son [14]: 1. Forma canónica controlable (First companion form) 2. Forma canónica observable (Second companion form) 3. Forma canónica de Jordan 4. Forma canónica controlable alternativa / Forma canónica Diagonal (Alternative first companion form, Toeplitz first companion form) En contraposición a los inconvenientes de la función de transferencia, el espacio de estados presenta una seria de ventajas [4], como se mencionó anteriormente: Aplicable a sistemas lineales y no lineales Permite analizar sistemas de más de una entrada y una salida Pueden ser variantes o invariantes en el tiempo Las condiciones iniciales no necesariamente deben ser nulas 14 2 Modelado en el Espacio de Estados Es capaz de proporcionar información de lo que está sucediendo en el interior del sistema en cada momento Los resultados que ofrece los presenta de una forma sencilla y elegante, lo cual es un punto a favor en cuanto a una mejor comprensión del estudio. Este punto parece de poca importancia, pero cuando se tratan sistemas complejos en los que se tienen múltiples datos interrelacionados entre sí, una buena visualización ayuda bastante a su mejor entendimiento. 2.4. No unicidad de la representación en espacio de estados de un sistema Como ya se viene advirtiendo en los apartados anteriores y se ha comprobado con los dos métodos analizados (Método de programación directa yMétodo de programación anidada), a un sistema descrito por su función de transferencia le corresponden al menos dos representaciones en espacio de estados diferentes (FCC yFCO respectivamente). Esto es debido principalmente al hecho de que la dinámica de un sistema puede ser descrita por multitud de variables de estado. Se pueden tomar variables de estado que sean combinaciones lineales de otras, que no alteren las propiedades del sistema, pero que sin embargo obliguen a modificar la estructura de resolución del problema. En este caso se obtendrían unas representaciones distintas pero equivalentes. Para demostrar este fenómeno, se utilizará una transformación mediante una matriz invertible P, la cual relacionará el actual vector de estado x(k)con otro ˜x(k)con variables de estado distintas mediante [13]: x(k) = P˜x(k)(2.24) Se obtendría una nueva ecuación de estado del sistema: P˜x(k+ 1) = GP ˜x(k) + Hu(k) =⇒˜x(k+ 1) = P−1GP ˜x(k) + P−1Hu(k)(2.25) Por lo que las ecuaciones que definen la nueva representación del sistema serían las siguientes: ˜x(k+ 1) = ˜ G˜x(k) + ˜ Hu(k) y(k) = ˜ C˜x(k) + ˜ Du(k)=⇒ ˜ G=P−1GP, ˜ H=P−1H ˜ C=CP, ˜ D=D(2.26) Se obtiene así un sistema en espacio de estados equivalente al anterior, pero con variables diferentes. 2.5. Resolución de las ecuaciones del espacio de estados Para la resolución de estas ecuaciones, es importante primero saber cuáles son los requisitos para su obtención. Éstas se pueden obtener mediante las ecuaciones diferenciales que representan un sistema. Los pasos a seguir a grandes rasgos son los siguientes [4] [17]: 1. Identificar las leyes o teorías que gobiernan el comportamiento que sigue el sistema (Leyes de termodinámica, Leyes dinámicas, Segunda Ley de Newton, Ley de voltajes y corrientes de Kirchoff, Ley de Ampere, Ley de Ohm, Ley de Boyle, etc.) 2. Seleccionar las variables de estado. Son las mínimas variables que determinan el comportamiento dinámico del sistema. En este paso se determina el orden ndel sistema que se va a estudiar. 15 2 Modelado en el Espacio de Estados 3. Encontrar la dinámica de cada estado. Es decir, se debe conocer cómo varía esa variable con respecto al tiempo, o lo que es lo mismo, su derivada con respecto al tiempo. En este paso son definidas las matrices A, B, C y D del sistema representado en 2.21. En este trabajo se explicarán tres métodos para obtener el valor del vector de estado. A partir del valor inicial x0, se obtendrá el valor para cualquier instante de tiempo k > 0, mediante alguno de los siguientes procesos: Procedimiento recursivo Matriz de transición de estados Método basado en la transformada Z 2.5.1. Procedimiento recursivo Este procedimiento debe su nombre al hecho de que es necesario realizar un proceso iterativo para obtener la solución. Si se realiza dicho proceso sobre las ecuaciones 2.3, las cuales representan un sistema LTI a partir de k= 0: x(1) = Gx(0) + Hu(0) x(2) = Gx(1) + Hu(1) = G2x(0) + GHu(0) + Hu(1) x(3) = Gx(2) + Hu(2) = G3x(0) + G2Hu(0) + GHu(1) + Hu(2) ... (2.27) Lo cual generalizado para cualquier k > 0quedaría de la siguiente forma: x(k) = Gkx(0) + k−1 X j=0 Gk−j−1Hu(j)(2.28) Simplemente observando la ecuación 2.28 se puede observar que los valores del actual x(k)dependerán tanto del estado inicial x0como de los valores de la entrada u(j). La expresión 2.29 determina la salida del sistema. y(k) = CGkx(0) + C k−1 X j=0 Gk−j−1Hu(j) + Du(k)(2.29) 2.5.2. Matriz de transición de estados Este proceso es un poco más limitado que el anterior, puesto que se presupone que no existe una señal de entrada u(k), y por tanto el estado actual solamente dependería del anterior. Convirtiendo esta información en una ecuación matemática se obtiene lo siguiente: x(k+ 1) = Gx(k)(2.30) Puesto que no posee una señal de entrada, se puede expresar la solución de la ecuación refiriéndose al estado inicial, puesto que éste es el único que se debe conocer para ir averiguando los siguientes estados actuales. Es decir, se necesita una función que, al conocer el momento en el que queremos calcular el estado, sea capaz de calcular dicho estado a partir del instante inicial. A esta función se le llamará ψ: 16 2 Modelado en el Espacio de Estados x(k) = ψkx(0) =⇒con :ψ(k+ 1) = Gψ(k)y ψ(0) = I=⇒es decir :ψ(k) = Gk(2.31) A dicha función ψkse le denomina matriz de transición de estados, y contiene toda la información sobre movimientos libres del sistema descrito por 2.30. Estos movimientos libres son de los que se hablaba anteriormente, y se refieren a los cambios de estado o su evolución en ausencia de la entrada uk. Las soluciones de la ecuación de estado y de la ecuación de salida son las siguientes: x(k) = ψ(k)x(0) + k−1 X j=0 ψ(k−j−1)Hu(j) =ψ(k)x(0) + k−1 X j=0 ψ(j)Hu(k−j−1) y(k) = Cψ(k)x(0) + C k−1 X j=0 ψ(j)Hu(k−j−1) + Du(k) (2.32) 2.5.3. Método basado en la transformada Z . Este método es complejo debido a que se necesita calcular la transformada de algunas de las expresiones. Como se ha dicho con anterioridad, la transformada Z difiere de la transformada de Laplace en que la primera se utiliza para sistemas que hacen uso de señales discretas (modelo espacio estados), mientras que el segundo se utiliza con señales continuas (función de transferencia, por ejemplo). Partiendo de 2.3 y realizando la transformada Z a ambos lados de la igualdad se obtiene: zX(z)−zx(0) = GX(z) + HU(z)⇔(zI −G)X(z) = zx(0) + HU(z) ⇔X(z)=(zI −G)−1zx(0) + (zI −G)−1HU(z)(2.33) cuya antitransformada sería la siguiente: x(k) = Z−1[(zI −G)−1z]x(0) + Z−1[(zI −G)−1HU(z)] (2.34) Si se compara esta expresión con la obtenida en el método recursivo 2.28, se pueden igualar algunos términos, con lo que quedaría: Gk=Z−1[(zI −G)−1z]y k−1 X j=0 Gk−j−1Hu(j) = Z−1[(zI −G)−1HU(z)] (2.35) 2.6. Linealización de las ecuaciones de estado La mayoría de los procesos que suceden en la naturaleza contienen un alto grado de no linealidad. Por este motivo, es importante que la ciencia y la técnica sean capaces de proporcionar métodos de resolución para este tipo de problemas. Puesto que la no linealidad de un fenómeno concreto tiene unas características muy particulares, desarrollar técnicas capaces de resolver un determinado tipo de 17 2 Modelado en el Espacio de Estados El concepto de observabilidad está relacionado con la condición de “observación” o estimación de las variables de estado a partir de las variables de salida, las cuales son generalmente medibles. Tiene que ver con la posibilidad de determinar el valor del vector de estados de un sistema a partir de observaciones de las salidas y las entradas de dicho sistema. “Se dice que un sistema es completamente observable si cada variable de estado del sistema afecta a alguna de las salidas”. [10] [13] Con el fin de aclarar estos términos, se ilustra el siguiente diagrama de bloques, y una explicación de la aplicación de la terminología sobre controlabilidad y observabilidad. Figura 2.5: a) Sistema de Control con realimentación del estado. b) Sistema de control con realimentación del estado y observador. La Figura 2.5 muestra un sistema con la dinámica descrita por la ecuación 2.55. Realimentando las variables de estado gracias a la matriz de ganancia K, descrita en la ecuación 2.56, se obtiene un sistema en bucle cerrado que se describe mediante 2.57. ˙x(t) = Ax(t) + Bu(t)(2.55) u(t) = −Kx(t) + r(t)(2.56) ˙x(t) = (A−BK)x(t) + Br(t)(2.57) Este procedimiento forman las bases del diseño por ubicación de polos mediante el proceso de realimentación. El objetivo sería encontrar la matriz Kde realimentación, tal que los valores característicos de (A−BK)tengan ciertos valores predeterminados a los que se pretende llegar. En este caso, se puede afirmar lo siguiente 1. El sistema de la ecuación 2.55 es controlable si existe una matriz de realimentación constante K que permite que los valores característicos de (A−BK)sean asignados de forma arbitraria. La controlabilidad juega un papel importante en la ubicación de polos en los sistemas de control. 2. Se puede dar el caso en el que no todas las variables de estado estén físicamente disponibles, y por tanto, se necesite implementar un “observador” que sea capaz de estimar el vector de estado a partir del vector de salida y(t), tal y como se muestra en la Figura 2.5. El vector ¯xse denomina vector de estado observado, y se usa para generar el control u(t)a través de la matriz de realimentación K. La condición de que tal observador pueda ser diseñado para el sistema se conoce como observabilidad del sistema. 24 2 Modelado en el Espacio de Estados La controlabilidad se puede definir tanto para los estados (Controlabilidad del estado) como para la salida (Controlabilidad de la salida). Existen diferencias en su definición, y se explica a continuación para un sistema de control en tiempo discreto (ya que se ha explicado como discretizar sistemas continuos en la sección anterior), lineal e invariante en el tiempo. 2.8.1. Controlabilidad del estado completo Se considera el sistema LTI de control definido por x((k+ 1)T) = Gx(kT) + Hu(kT)(2.58) donde x(kT) = vector estado (dimensión n) en el k-ésimo instante de muestreo u(kT) = señal de control en el k-ésimo instante de muestreo. Se supone constante para kT ≤t < (k+T)T G= matriz de n×n H= matriz de n×1 T = período de muestreo Se dice que un sistema es de estado completamente controlable, o simplemente de estado controlable, si existe una señal de control constante por intervalos u(kT)definida a lo largo de un número finito de períodos de muestreo de forma que, partiendo de un estado inicial, el estado x(kT) puede ser transferido al estado deseado xfen nperíodos de muestreo como máximo. Se deduce a continuación la condición para la controlabilidad completa del estado x(nT) = Gnx(0) + n−1 X j=0 Gn−j−1Hu(jT) =Gnx(0) + Gn−1Hu(0) + Gn−2Hu(T) + ... +Hu((n−1)T) =⇒ x(nT)−Gnx(0) = [H. . .GH. . .· · · . . .Gn−1H]     u((n−1)T) u((n−2)T) . . . u(0)      (2.59) donde la matriz Mc= [H. . .GH. . .· · · . . .Gn−1H]se denomina matriz de controlabilidad. Puesto que Hes una matriz n×1, se tiene que cada una de las matrices que componen la matriz de controlabilidad es una matriz de n×1o un vector columna. La condición necesaria y suficiente para que el sistema sea completamente controlable es Rango[H. . .GH. . .· · · . . .Gn−1H] = n(2.60) Se puede demostrar que esta condición también es válida para un sistema en el que u(kT)sea un vector de dimensión r. [13] 25 2 Modelado en el Espacio de Estados 2.8.2. Controlabilidad de la salida Otro caso posible que se puede encontrar a la hora de diseñar un sistema de control, es preferir controlar la salida en vez del estado. Este tipo de controlabilidad no necesita que se produzca una controlabilidad completa del estado, y por esta razón, es necesaria definirla por separado. Es uno de los objetivos más comunes, controlar la evolución de la salida del sistema. En base a las ecuaciones x((k+ 1)T) = Gx(kT) + Hu(kT) y(kT) = Cx(kT) + Du(kT)(2.61) donde x(kT)= vector estado (dimensión n) en el k-ésimo instante de muestreo. u(kT)= señal de control (dimensión r) en el k-ésimo instante de muestreo y(kT)= vector de salida (dimensión m) en el k-ésimo instante de muestreo G= matriz de n×n H= matriz de n×r C= matriz de m×n D= matriz de m×r Se dice que el sistema 2.61 es de salida completamente controlable (o simplemente salida controlable) si es posible tener una señal de control no restringida u(kT), definida en un conjunto finito de períodos de muestreo 0≤kT ≤nT tales que, partiendo de la salida inicial y(0), la salida y(kT)pueda ser transferida al punto deseado (punto arbitrario) yfen el espacio de salidas, en n períodos de muestreo como máximo. Se deduce a continuación las condiciones necesarias para que se produzca la controlabilidad de la salida. La salida puede ser dada por y(nT) = Cx(nT) + Du(nT) =CGnx(0) + n−1 X j=0 CGn−j−1Hu(jT) + Du(nT)(2.62) obteniéndose y(nT)−CGnx(0) = n−1 X j=0 CGn−j−1Hu(jT) + Du(nT) =CGn−1Hu(0) + CGn−2Hu(T) + ... +CHu((n−1)T) + Du(nT) = [D. . .CH. . .CGH. . .· · · . . .CGn−1H]    u(nT) u((n−1)T) ... u(0)     (2.63) Una condición necesaria y suficiente para que el sistema sea de salida completamente controlable es que la matriz de la última expresión de 2.63 sea de rango m, es decir 26 2 Modelado en el Espacio de Estados Rango[D. . .CH. . .CGH. . .· · · . . .CGn−1H] = m(2.64) Cabe destacar que la existencia de la matriz D en la ecuación de salida ayuda a establecer la controlabilidad del sistema. En el caso de que no existiera la matriz D (salida y(kT) = Cx(kT)), la condición necesaria y suficiente sería Rango[CH. . .CGH. . .· · · . . .CGn−1H] = m(2.65) En este caso, la controlabilidad de la salida se produce si y sólo si los m renglones de C son linealmente independientes. Si se realiza un análisis de ambos casos estudiados, la forma en la que aparece Dtiene más probabilidades de ser controlable. Al introducirse una columna extra en la matriz de controlabilidad (la correspondiente a D), se podría dar la situación en la que se pase de tener m−1columnas linealmente independientes a m, por lo que se conseguiría la controlabilidad. Por simplificar, encontrar mvectores linealmente independientes es igual o más fácil entre n+ 1 vectores que en solo nvectores. [13] 2.8.3. Observabilidad Para analizar la observabilidad de un sistema de control, se hará uso de un sistema LTI discreto descrito por x((k+ 1)T) = Gx(kT) y(kT) = Cx(kT )(2.66) donde x(kT) = vector de estado (dimensión n) en el k-ésimo instante de muestreo y(kT) = vector de salida (dimensión m) en el k-ésimo instante de muestreo G= matriz de n×n C= matriz de m×n El hecho de utilizar en este caso un sistema sin excitación es que, para la investigación de la condición necesaria y suficiente para la completa observabilidad, basta considerar el sistema 2.66, puesto que las matrices que se añaden e un sistema con variables de entrada son conocidas y solamente supondrían “complejidad” al estudio de las ecuaciones. Se dice que un sistema es completamente observable si cualquier estado inicial x(0) puede determinarse a partir de la “observación” de y(kT)sobre un número finito de períodos de muestreo, es decir, si cualquier transición del estado de manera eventual afecta a todos los elementos del vector de salida. Otra posible definición sería: un sistema es completamente observable si cada variable de estado del sistema afecta alguna de las salidas. La observabilidad es un fenómeno realmente útil para dar solución a un problema en el que las variables de estado de un sistema no son medibles puesto que dichas variables no son accesibles. La solución de la ecuación 2.66 será x(kT) = Gkx(0) =⇒por tanto y(kT) = CGkx(0) (2.67) 27 2 Modelado en el Espacio de Estados La completa observabilidad implica lo siguiente: conocidos y(0), y(T), y(T), ..., es posible determinar x1(0), x2(0), ..., xn(0). Puesto que se necesitan encontrar nincógnitas, se necesitan únicamente nvalores diferentes de y(kT ). Esto significa que se pueden utilizar los primeros nvalores de y(kT)ó y(0), y(T), y(T), ..., y((n−1)T)que permiten determinar x1(0), x2(0), ..., xn(0). Es decir, dados y(0) = Cx(0) y(T) = CGx(0) . . . y((n−1)T) = CGn−1x(0) se debe ser capaz de determinar x1(0), x2(0), ..., xn(0). Realizando el análisis dimensional correspondiente, se puede observar que y(kT )es un vector de dimensión m, y puesto que se tienen necuaciones, se llega a un sistema de n×mecuaciones, todas ellas incluyendo x1(0), x2(0), ..., xn(0). Si se fija como objetivo obtener una solución de x1(0), x2(0), ..., xn(0) a partir de esas n×mecuaciones, se debe imponer que exactamente nde ellas sean linealmente independientes, o lo que es lo mismo Rango[C∗. . .G∗C∗. . .· · · . . .(G∗)n−1C∗] = n(2.68) donde C∗significa la transpuesta conjugada de C. Se ha llegado, por tanto, a la condición necesaria y suficiente para que el sistema de las ecuaciones 2.66 sea completamente observable. La matriz de 2.68 se conoce como matriz de observabilidad. [13] 2.8.4. Principio de Dualidad Ahora que se han definido los conceptos de controlabilidad y observabilidad, se analizará a continuación la relación que existe entre ambas definiciones. Para ello, se considerarán los dos sistemas siguientes x((k+ 1)T) = Gx(kT) + Hu(kT) y(kT) = Cx(kT)(2.69) el cual será nombrado como S1. Para S2se tendrá lo siguiente ˆx((k+ 1)T) = G∗ˆx(kT ) + C∗ˆu(kT) ˆy(kT) = H∗ˆx(kT)(2.70) donde x(kT) y ˆ x(kT) = vector de estado (dimensión n) en el k-ésimo instante de muestreo u(kT) y ˆ u(kT) = vector de control (dimensión r) en el k-ésimo instante de muestreo y(kT) y ˆ y(kT) = vector de salida (dimensión m) en el k-ésimo instante de muestreo G= matriz de n×n,G∗= transpuesta conjugada de G H= matriz de n×r,H∗= transpuesta conjugada de H 28 2 Modelado en el Espacio de Estados C= matriz de m×n,C∗= transpuesta conjugada de C La analogía que existe entre controlabilidad y observabilidad se denomina Principio de dualidad, y se conoce gracias a Kalman. Y dice así “El sistema S1definido por las ecuaciones 2.69 es de estado completamente controlable (observable), si y sólo si el sistema S2definido por las ecuaciones 2.70 es de estado completamente observable (controlable)”. La demostración de este principio es bastante sencilla y se basa en las condiciones necesarias y suficientes para la controlabilidad y observabilidad completas. Si se escriben ambas condiciones, analizando ambos sistemas por separado, se comprenderá mejor 1. Controlabilidad S1: Rango[H. . .GH. . .· · · . . .Gn−1H]=n 2. Observabilidad S1: Rango[C∗. . .G∗C∗. . .· · · . . .(G∗)n−1C∗]=n 3. Controlabilidad S2: Rango[C∗. . .G∗C∗. . .· · · . . .(G∗)n−1C∗]=n 4. Observabilidad S2: Rango[H. . .GH. . .· · · . . .Gn−1H]=n Como se puede observar comparando ambas condiciones, se evidencia la verdad del principio de dualidad. Mediante la utilización de este principio, la observabilidad de un sistema dado puede verificarse al probar la controlabilidad del estado de su dual. [13] 2.9. Transformación de un sistema en formas canónicas Existen una serie de las llamadas transformaciones de similitud las cuales son realmente útiles a la hora de trabajar en el análisis y el diseño en el dominio de estado. Es posible que se quiera trabajar con este tipo de ecuaciones particulares por diversas razones, como se explicarán a continuación. Se denominan transformaciones de similitud puesto que tanto el sistema de partida como el sistema transformado conservan las mismas ecuaciones características, vectores característicos, valores característicos y la función de transferencia. A continuación se describen las transformaciones de la Forma Canónica Controlable (FCC) y de la Forma Canónica Observable (FCO) a partir del siguiente sistema x(k+ 1) = Gx(k) + Hu(k) y(k) = Cx(k) + Du(k)(2.71) 2.9.1. Forma Canónica Controlable Considerese la matriz de transformación P=SM, donde S=hH. . .GH . . .G2H. . .· · · . . .Gn−1Hi, M =        a1a2· · · an−11 a2a3· · · 1 0 . . .. . ..... . .. . . an−11· · · 0 0 1 0 · · · 0 0        (2.72) 29 2 Modelado en el Espacio de Estados donde Ses conocida como la matriz de controlabilidad yaison los coeficientes de la ecuación característica de G, la cual es |sI −A|=sn+an−1sn−1+· · · +a1s+a0= 0 (2.73) Definiendo el estado x(k)a partir de la matriz de transformación Pen función de otro estado ˆx(k) x(k) = Pˆx(k)(2.74) Se puede afirmar que el sistema ˆx(k+ 1) = ˆ Gˆx(k) + ˆ Hu(k) y(k) = ˆ Cx +ˆ Du(k)(2.75) se encuentra en Forma Canónica Controlable (FCC) si ˆ G=P−1GP =        010· · · 0 001· · · 0 . . .. . .. . ..... . . 000· · · 1 −a0−a1−a2· · · −an−1        ,ˆ H=P−1H=        0 0 . . . 0 1        ˆ C=CP, ˆ D=D (2.76) donde las matrices ˆ Cyˆ Dno siguen ningún patrón en particular. Este tipo de transformaciones tiene un requisito indispensable: que la matriz P−1exista, lo cual implica que la matriz Stenga inversa, puesto que la inversa de Msiempre existe (su determinante es (−1)n−1, nunca cero). [10] [13] Esta forma de expresar en modo de espacio estados una función de transferencia garantiza que el sistema que se esté modelando sea controlable, término que se ha descrito con anterioridad. Cuando un sistema es controlable, es decir, se pueden modificar todos y cada uno de los estados mediante las entradas, se puede expresar mediante su FCC. En el caso del estudio de la trayectoria de un vehículo no tripulado, interesa que el modelo que simula el sistema sea controlable, y por tanto, se podría expresar según su FCC. Esa controlabilidad será la que permita al piloto modificar la velocidad, posición, actitud y cualquier otra variable definida como estado del sistema. 2.9.2. Forma Canónica Observable Este método es bastante similar, y se le considera una forma dual de la transformación de la FCC. El sistema de 2.71 se transforma a la Forma Canónica Observable mediante la matriz Q x(k) = Qˆx(k), Q = (MV )−1(2.77) donde Mviene definido por la ecuación 2.72, y V=        C CG CG2 . . . CGn−1        (2.78) 30 2 Modelado en el Espacio de Estados V se denomina comúnmente matriz de observabilidad, y debe cumplir el requisito de que V−1 debe existir para que la transformación FCO sea posible. Las distintas matrices involucradas en las ecuaciones transformadas quedarían de la siguiente forma ˆ G=Q−1GQ =        0 0 · · · 0−a0 1 0 · · · 0−a1 0 1 · · · 0−a2 . . .. . ..... . .. . . 0 0 · · · 1−an−1        ,ˆ C=CQ =0 0 · · · 0 1  ˆ H=Q−1H, ˆ D=D (2.79) donde los elementos de las matrices ˆ Hyˆ Dno están restringidos de ninguna forma, sino que dependerán de la función de transferencia de partida. Cabe destacar que si se analizan estas matrices de ˆ Gyˆ Cse puede observar que equivalen a las traspuestas de las matrices ˆ Gyˆ Hde la expresión 2.76, y por ello se le denominaba forma dual. [10] [13] Para poder expresar un modelo mediante su FCO es necesario que el sistema bajo estudio sea completamente observable, es decir, que no haya salidas que no dependan de los estados directamente, o dicho de otra forma, que cada variable de estado afecte a alguna salida. La principal ventaja que poseen estas dos representaciones, FCC y FCO, son principalmente la rapidez con la que pueden ser calculadas. Una vez que se conozca que el sistema modelado cumple con los requisitos de controlabilidad u observabilidad, obtener el modelo de espacio estados a través de su ecuación diferencial o función de transferencia es inmediato, simplemente mirando los coeficientes. En el caso de no cumplir con dichos requisitos, sería necesario realizar cálculos más tediosos para obtener las matrices del modelo de espacio estados. 2.10. Descripción de un sistema en parte controlable/observable y no controlable/no observable En la práctica, es posible que el problema bajo estudio no sea completamente controlable/observable, es decir, que tengan elementos del sistema que sí lo son pero otros no. Se verá en este apartado que es posible descomponer un sistema en parte controlable/no controlable y observable/no observable. Esta idea resulta de utilidad para un mejor análisis de los datos de un sistema, pudiendo resolver la parte controlable/observable por un método y estudiar minuciosamente la parte restante del sistema. Se considerarán las mismas ecuaciones de estado que se han venido estudiando hasta el momento, como las de la expresión 2.71. 2.10.1. Parte controlable/no controlable Una condición que se presupone para que en un sistema exista parte no controlable, es que su rango sea menor que la dimensión de la ecuación de estado, como se demostró en las condiciones necesarias y suficientes para la controlabilidad. Es decir Rango H GH · · · Gn−1H=n1< n (2.80) De nuevo se usará la matriz Pcomo matriz de transformación o cambio de coordenadas, de dimensión n, definida como P−1=p1p2· · · pn1· · · pn(2.81) 31 2 Modelado en el Espacio de Estados donde las n1primeras columnas son linealmente independientes para que se cumpla la condición de la expresión 2.80. El resto se eligen arbitrariamente con el objetivo de hacer Puna matriz no singular. Esta transformación (¯x=Px) lleva al sistema a la siguiente expresión ˙ ¯xC ˙ ¯x¯ C=¯ GC¯ G12 0¯ G¯ C xC x¯ C+¯ HC 0u y=¯ CC¯ C¯ Cx+Du (2.82) Los estados en las nuevas coordenadas se descomponen en ¯xC:n1estados controlables ¯x¯ C:n−n1estados no controlables Figura 2.6: Descomposición parte controlable y no controlable La ecuación de estado de orden reducido (puesto que se le eliminan aquellos elementos no controlables) de los estados controlables quedaría ˙ ¯xC=¯ GC¯xC+¯ HCu ¯y=¯ CC¯x+Du (2.83) el cual es controlable y tiene la misma función de transferencia que el estado original. Como dato interesante, la función de Matlab ctrbf transforma una ecuación de estado en su forma canónica controlable/no controlable. 2.10.2. Parte observable/no observable Este fenómeno es bastante parecido al anterior, pero tiene algunas diferencias debido a la dualidad mencionada previamente. En este caso, se supondrá que el rango será Rango      C CG . . . CGn−1      =n2< n, P =           p1 p2 . . . pn2 . . . pn           (2.84) donde Pvuelve a ser la matriz para cambiar de coordenadas y sus primeras n2filas son linealmente independientes, y el resto se eligen arbitrariamente para que no sea singular. De nuevo, con la transformación ¯x=Px el sistema quedaría 32 2 Modelado en el Espacio de Estados ˙ ¯xO ˙ ¯x¯ O=¯ GO0 ¯ G21 ¯ G¯ O ¯xO x¯ O+¯ HO H¯ Ou y=¯ CO0¯x+Du (2.85) Los estados en las nuevas coordenadas se descomponen en ¯xO:n2estados observables ¯x¯ O:n−n2estados no observables Figura 2.7: Descomposición parte observable y no observable La ecuación de estado de orden reducido (puesto que se le eliminan aquellos elementos no observables) de los estados observables quedaría ˙ ¯xO=¯ GO¯xO+¯ HOu ¯y=¯ CO¯x+Du (2.86) el cual es observable y tiene la misma función de transferencia que el estado original. Como dato interesante, la función de Matlab obsvf transforma una ecuación de estado en su forma canónica observable/no observable. 33 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Rentabilidad y prototipado. Requiere menos hardware en comparación con los prototipos físicos, de modo que supone menor coste y son más rápidos de construir. Fidelidad y verosimilitud. Alcanzan niveles de fidelidad inalcanzables por simulaciones puramente virtuales, puesto que estos no son capaces de reproducir completamente algunas de las características o atributos del sistema real. Velocidad de simulación. Alcanzan mayor velocidad que las simulaciones virtuales de los mismos fenómenos (por ejemplo, simulaciones de motor IC basados en Dinámica de Fluidos Computacional). Repetitividad. Aquellos sistemas que normalmente operan en entornos muy variables (por ejemplo, los sistemas de suspensión de vehículos fuera de la carretera) a menudo pueden ser probados en entornos de laboratorio controlados a través de la simulación HIL, que pueden aumentar significativamente la repetibilidad. No altera la naturaleza. Hace posible la simulación de eventos destructivos (accidentes de vehículos, interceptación de misiles, etc.) sin incurrir en una destrucción real, y por lo tanto, costosa. Integralidad. Hace posible la simulación de un sistema en un rango mucho más amplio de sus condiciones de funcionamiento posibles, a través de la creación de prototipos puramente físicos. Seguridad. Se pueden utilizar para entrenar a los operadores humanos (por ejemplo, pilotos de avión) o a sistemas críticos para la seguridad (por ejemplo, aviones supersónicos) en entornos significativamente más seguras (por ejemplo, simuladores de vuelo). Coexistencia de sistemas de ingeniería. Permite que diferentes equipos desarrollen las diferentes partes de un sistema en el hardware sin perder de vista los problemas de integración, permitiendo de este modo la ingeniería de sistemas concurrentes. Por otra parte, el sistema SIL corren también en tiempo real con RTOS (del inglés, Real Time Operative Systems) y sirven para desarrollar leyes de control y de navegación, donde todo se simula por software (tanto los sensores como los actuadores). En el caso de los vehículos aéreos no tripulados, el simulador SIL podría consistir en varios ordenadores en los que se simula la GCS (del inglés, Ground Control Station), y un ordenador en el que se simula el modelo de aeronave (modelo dinámico, aerodinámico, de masa e inercial, flight control: navegación-guiado-control, el sistema de gestión de la misión, generación y envío de las medidas tomadas por los sensores y actuadores, etc.). En cualquier caso, ambas técnicas son utilizadas para llevar a cabo experimentos que sean lo más similares posibles a las actuaciones reales que tendrán que realizar los sistemas embarcados. De esta forma, se pueden analizar los fallos o carencias que posean los dispositivos antes de ponerlos en funcionamiento, por lo que supondrá un ahorro considerable en el diseño de sistemas. 3.3. Sensores para medir distancias y proximidad Se describen este tipo de sensores en primer lugar debido a su importancia. Dentro de las distintas funcionalidades que puede tener un UAV, en la mayoría de ellas será necesario conocer los objetos que lo rodean para realizar una función concreta. Por ejemplo, en campañas anti-incendios se necesita conocer el área invadida por las llamas, en control de plagas en agricultura es necesario cubrir todo el terreno con los correspondientes agentes plaguicidas, en mensajería aérea también es preciso conocer el punto exacto de entrega, etc. Estos sensores de proximidad surgen de esta necesidad, de conocer la posición en la que no existe contacto entre el actuador y el detector. Para todas aquellas funciones, existen los sensores que miden la proximidad y la presencia de objetos situados a una distancia máxima de alcance del sensor. Se describirán en primer lugar sensores 40 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados que hacen uso de características eléctricas (capacitivos); después se explicarán aquellos que usan propiedades eléctricas y magnéticas (inductivos, basados en efecto Hall); y por último los que hacen uso de propiedades ópticas y ultrasonidos. 3.3.1. Sensores capacitivos Cabe destacar en primer lugar la diferencia existente entre sensores capacitivos e inductivos. Los capacitivos funcionan detectando los cambios que se producen en la capacidad parásita que se origina entre el detector y objeto que se desea medir. Son utilizados para medir la distancia a objetos metálicos y no metálicos, como pueden ser la madera, líquidos y plásticos. Por su parte, los sensores inductivos están basados en el amortiguamiento que se produce en el campo magnético a causa de corrientes inducidas (o corrientes de Foucault) de los materiales situados alrededor del sensor. En este caso el material debe ser metálico para su funcionamiento. Los sensores capacitivos se basan en el esquema representado en la Figura 3.4. Se observan condensadores genéricos, en los cuales se dispone de dos placas metálicas llamadas armaduras, separadas por un dieléctrico. Figura 3.4: Esquemas de sensores capacitivos [22] Los elementos capacitivos utilizados como condensadores son variables, de forma que el desplazamiento a medir provocará un desplazamiento en algún componente del condensador, y por lo tanto, una modificación de su capacidad. En la Figura 3.4 se observa cómo al acercar un objeto se desplaza un componente en el condensador y por lo tanto, se modifica la capacidad. Esa variación permite medir el desplazamiento sufrido o la cercanía con el objeto. Entre las dos placas se almacena una electricidad, la cual puede verse variada por la modificación de la posición del dieléctrico o la disposición de las mismas. En el caso de placas paralelas entre sí, el campo eléctrico creado entre ellas es uniforme, y la diferencia de potencial es VR=1 a Qd S(3.2) donde Qes la carga de cada lámina, Sel área, ala constante dieléctrica del medio, y dla separación. La capacidad del condensador es la siguiente C=a S d(3.3) El cálculo de la distancia a medir se basa en el principio de que el potencial que se almacena es inversamente proporcional a la distancia que separa las placas del condensador. La cercanía con el objeto provocará un aumento de la capacitancia. Dicho aumento, mediante el proceso de calibrado conveniente, es traducido a una distancia como señal de salida del sensor. Para medir distancias mayores a algunos milímetros, la sensibilidad de estos sensores disminuye notablemente. En la Figura 3.5 se puede observar cómo la cercanía de un objeto, sea o no conductivo, implica un aumento del campo eléctrico. 41 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Figura 3.5: Esquema de sensor capacitivo como detector de proximidad Una forma fácil de simplificar este funcionamiento sería: un cambio en la distancia se traduce en un desplazamiento en algún componente que provoca cambio de capacidad, mediante la variación del campo eléctrico, y esa capacidad es traducida a un cambio en la salida del sensor para su posterior interpretación. 3.3.2. Sensores inductivos Para utilizar este tipo de sensores es necesario que el objeto del cual se quiere medir su distancia sea ferromagnético. Los sensores poseen una bobina y un imán permanente. Al colocar el iman cerca de la bobina y cambiar de posición el objeto ferromagnético, entre ambos componentes se producirá una variación del flujo magnético a través de la bobina, lo que inducirá una fuerza electromotriz. Esta corriente es la que se utilizará como detector de presencia del objeto a medir. Normalmente son utilizados cuando el desplazamiento relativo entre el sensor y el objeto es lineal. Un ejemplo de sensores que utilizan este principio son los llamados transformadores diferenciales. Se dispone de una bobina primaria central y dos bobinas secundarias laterales enrolladas sobre un núcleo magnético, de forma que los desplazamientos que van a medir estos sensores van a ser los de este núcleo. Se proporciona una intensidad y tensión conocidas a la primaria, la cual es transformada en las bobinas secundarias en intensidades y tensiones diferentes. Figura 3.6: Ejemplo de transformador diferencial [22] En la Figura 3.6 se puede observar el funcionamiento del transformador diferencial. El desplazamiento producido en el núcleo magnético se produce de tal forma que alguna de las bobinas secundarias no cubra completamente el núcleo. Como consecuencia de esto, la corriente inducida en una secundaria será mayor que en la otra, y a partir de esta diferencia, se podrá medir el desplazamiento sufrido por el núcleo. Con el núcleo centrado, la salida será de 0V, mientras que si se desplaza a algún lado habrá más o menos corriente inducida en una de las dos secundarias, la cual se traducirá en distancia de desplazamiento. 42 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Figura 3.7: Sensores capacitivos e inductivos reales (cortesía de Siriored) Este tipo de sensores, a pesar de no poder medir grandes desplazamientos, tienen ventajas como su rápida respuesta, alta resolución, linealidad, baja histéresis y repetitividad. 3.3.3. Basados en efecto Hall Para utilizar un sensor que se base en este principio, es necesario que el objeto a medir sea capaz de producir un campo magnético. Se obtendrá una referencia de su cercanía o del desplazamiento que ha sufrido a partir de la diferencia de potencial generada como consecuencia del efecto Hall. Esta diferencia de potencial será mayor cuanto más intenso sea el campo magnético o más próximo se encuentre el objeto a medir. ¿Cómo funcionan realmente este tipo de sensores? Utilizando la presencia de un campo magnético en un semiconductor para producir cambios en la corriente eléctrica generada. Se utiliza esta presencia o ausencia del campo magnético para proporcionar un determinado nivel de tensión V, como se muestra en la Figura 3.8. Se observa que el campo magnético es perpendicular a la corriente eléctrica que atraviesa la placa conductora. Se genera así el campo eléctrico debido a la polarización de la placa en lado positivo y negativo, compensando el campo magnético. Figura 3.8: Efecto Hall [22] Entre ambos extremos de la placa se genera la siguiente diferencia de potencial V=KH Bfi d(3.4) donde KHes el coeficiente de Hall, Bfes la densidad del flujo magnético, ies la intensidad de corriente y del grosor de la placa. Típicamente este tipo de sensores se encuentran instalados en semiconductores. La electrónica integrada que poseen les permite proporcionar una señal que se encuentra amplificada y condicionada, lo cual supone una actuación más directa, fácil y económica. En [5], proyecto en el que se intenta mejorar el sistema anemométrico y de monitorización energética en un UAV (Proyecto Céfiro del Departamento de Aeroespacial de la Universidad de Sevilla), se utilizan este tipo de sensores basados en el efecto Hall. En dicho proyecto, se contaban los pasos de las bobinas de un motor eléctrico por un punto estático mediante un sensor de efecto Hall, para la medición del régimen de giro del motor. En este caso se utilizó un Allegro 1120 EUA-T. Dicho sensor actuaba como interruptor, el cual se accionaba debido al paso del campo magnético generado por cada bobina del motor, puesto que esta elevaba la tensión terminal del sensor de 0Va 5V. De esta forma era capaz de contar el tiempo entre vueltas, y así, la velocidad de giro del motor. 43 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados 3.3.4. Basados en ultrasonidos Los ultrasonidos son ondas con la capacidad de ser reflejadas si encuentran en su camino una discontinuidad o algún elemento extraño en el medio en el que viajan. Son prácticamente iguales a las ondas sonoras, por lo que viajan a una velocidad de 340m/s, pero poseen una mayor frecuencia (en torno a los 20kHz, por encima del umbral del oído humano, 16kHz). Su funcionamiento es “sencillo”: el emisor lanza un tren de pulsos ultrasónicos y espera el rebote en el receptor, para medir el tiempo transcurrido para así calcular la distancia a la que se encuentra. La reflexión de la onda en el objeto a medir se produce debido a la diferencia de impedancias acústicas entre el medio y el objeto. El rango de actuación, en cuanto a distancia que pueden medir estos sensores, es mayor del que poseen los presentados con anterioridad. Pero también cabe destacar que presentan una serie de problemas de implementación bastante comunes, los cuales dificultan las buenas mediciones de distancias lejanas. Algunos de los problemas más comunes son: Ángulo de incidencia. La dirección del reflejo de la onda depende directamente de este ángulo. Cuanto menor sea, más probable será que no se detecte el eco o que se produzcan medidas erróneas. Figura 3.9: Explicación del error producido debido al ángulo de incidencia Superficie. Una superficie lisa agrava el problema del ángulo de incidencia. Cuanta más rugosidad, más superficie donde la onda podrá rebotar y por tanto, más probabilidad de que la lectura sea correcta. Se suele concluir que para que se produzca una buena reflexión, las irregularidades deben ser del orden de magnitud de la longitud de onda del ultrasonido. Ambiente. Las turbulencias debidas a las corrientes de aire pueden dificultar la detección del ultrasonido. La temperatura provoca cambios en la densidad del aire, y por tanto, en la velocidad con la que se propaga la onda, con el correspondiente error en la medición de la distancia. Cercanía. Algunos de estos sensores necesitan un intervalo de tiempo desde que envían la señal para preparar al receptor. Si el eco rebota antes de ese tiempo, no leerá la señal. Se suele decir que se necesita una distancia mínima para la detección de distancias debido a este tiempo mínimo de espera. Rango de detección. El campo de acción de la onda ultrasónica tiene una forma cónica, de forma que solamente se pueden obtener datos de los objetos dentro de dicho cono. Este hecho limita el rango de detección del sensor y además, supone una incertidumbre debido a que no se especifica la localización angular dentro de dicho cono. Crosstalk. Detección de falsos ecos, producida en aquellas áreas en las que se usan distintos sensores de ultrasonidos. Se produce crosstalk cuando una señal es recibida por un sensor distinto del emisor. También se producen errores cuando la onda sufre varias reflexiones antes de volver al receptor, lo cual proporciona una medida errónea de la distancia. 44 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Figura 3.10: Sensores de efecto Hall (izquierda, Bernstein) y ultrasonidos (derecha, Arduino HCSR04 [9]) 3.3.5. Sensores de Espectro Infrarrojo Este tipo de sensores a menudo es definido en apartados de sensores de luz, pero su principal aplicación reside en el seguimiento o evitación de obstáculos, por lo que se ha decidido incluirlo en este apartado. Los sensores infrarrojos utilizan parte del espectro denominado infrarrojo, invisible para el ser humano. A diferencia de otros sensores de luz, estos presentan una baja proporción de interferencias. No proporcionan información directa de la distancia hacia un objeto en general (existen modelos como GP2D02 y GP2D12 de la empresa Sharp que sí lo hacen), pero sí indican información acerca de si hay obstáculo o no en su cono de detección, que suele ser más estrecho que el de ultrasonidos. Constan de una fuente luminosa (lámparas, diodos LED, láser, etc.) y una célula encargada de la recepción de la señal, que puede ser un fotodiodo o un fototransistor (ambos explicados en el apartado de sensores de luz). Suelen ser útiles para distancias de decenas de milímetros, o incluso algunos centímetros. El sensor de la Figura 3.11 de la empresa Sharp consigue medir hasta unos 24cm, lo cual no lo convierte en el más adecuado para su uso en UAVs, puesto que se trata de una distancia de muy corto alcance. Es necesario que los sistemas embarcados sean capaces de detectar objetos en largas distancias, para disponer del tiempo necesario para tomar una decisión y llevarla a cabo. Es posible que se use como sistema de visión de forma complementaria, como por ejemplo para actuaciones en el interior de edificios para rastreado de grietas en muros, pero siempre suelen estar integrados con algún otro de más alcance. Este es el fundamento del fenómeno “Sense and Avoid”, explicado más en detalle en el siguiente apartado de visión artificial. Figura 3.11: Sensores de infrarrojos, Sharp GP2D15) Algunas de las aplicaciones donde su uso es común es en el seguimiento de un trazado de líneas, de paredes y detección de obstáculos en un circuito de baja velocidad. Tiene el inconveniente de ser sensible a la luz ambiente y a la reflectividad de los objetos. 3.3.6. Visión artificial Cada vez más se pretende que los vehículos aéreos no tripulados sustituyan a los seres humanos en distintas labores de localización y reconocimiento, ya sea por medidas de seguridad, rapidez, eficacia o una combinación de todas ellas. Para que este fenómeno sea posible, es necesario dotar de una 45 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados visión tridimensional al UAV, puesto que el entorno en el que trabajarán tendrá dichas características espaciales. Aparte de esto, también será necesario implementar algún tipo de sistema que sea capaz de adaptarse a los cambios en sus alrededores, puesto que el entorno de trabajo real, fuera de laboratorio, es dinámico. A su vez, el ámbito en el que se trabaja suele ser de tipo no estructurado, es decir, no se tienen conocimientos a priori y precisos de cómo se van a comportar los elementos que lo conforman. Por ejemplo, en el caso de reconocimiento de cultivos, se tienen conocimientos de lo que se espera ver, pero no se trata de una planta industrial en la que se conoce cómo van a actuar los robots manipuladores en cada momento. Es por este motivo que es necesario dotar a los drones con capacidades sensoriales lo más similares posibles a las del ser humano dentro de lo posible. Para aclarar el concepto, [22] hace un símil entre los elementos del ojo biológico y el sistema visual artificial: Sistema de entrada de información equivalente al ojo: estaría formado por un conjunto de cámaras de vídeo y de algún sistema de adquisición de imágenes. Más adelante serán procesadas. Sistema de almacenamiento y procesamiento equivalente al cerebro: formado por una tarjeta de procesamiento y algún ordenador cuya CPU se encargue de realizar aquellos procesos que la tarjeta no esté capacitada. Sistema de salida/visualización (sin equivalente en el ojo biológico): suele ser un monitor para mostrar la información que ha sido procesada. Lo más recomendable es que sea en tiempo real para una mayor rapidez de actuación en caso de algún imprevisto. Existe una amplia variedad de cámaras que pueden ser implementadas en UAV. Dependiendo de la función que vaya a realizar, de la capacidad de soportar carga del vehículo, de la resolución que se necesite, del precio que se esté dispuesto a pagar ..., se podrá escoger la más adecuada. Existe un gran interés en desarrollar software y hardware capaces de proporcionar sistemas fiables dentro del ámbito del llamado “Sense and Avoid”. Este concepto hace referencia a la integración de todo tipo de sistemas capaces de realizar tareas similares a las que puede realizar un ser humano en cuanto a temas de seguridad y vigilancia. Para ello, es necesario el uso de sensores de visión artificial capaces de detectar obstáculos y calcular a qué distancia se encuentran para así evitar colisiones indeseadas. Se pretende mejorar la calidad del “Sense and Avoid” primero con vehículos aéreos no tripulados de pequeño tamaño, pero con el objetivo de en un futuro implementarlo en grandes aeronaves para reforzar el trabajo de los pilotos o incluso, en un futuro no se sabe cómo de lejano, sustituirlos. Existen dos funciones claramente diferenciadas: Sense. Hace referencia a la observación del llamado “intruso”, para obtener la mayor información posible acerca de sus características y régimen de vuelo. Como objetivos prioritarios se encuentran conocer el posicionamiento, rumbo y velocidad. En vuelos tripulados visuales VFR (del inglés, Visual Flight Rules), esta función se realiza a simple vista con ayuda de una simple radio. En vuelos instrumentales IFR (del inglés, Instrumental Flight Rules) o nocturno VFRN (del inglés, Visual Flight Rules Night) se precisan sistemas de ayuda a la visión del piloto, como son estos sensores de visión artificial. Avoid. Esta función es la encargada de analizar y procesar la información recibida del Sense. Debe tener la capacidad de decidir si el intruso detectado es conflictivo o no, si existe riesgo de colisión o no, y en caso afirmativo advertir al piloto. En el caso de los vehículos no tripulados, debe enviar información a la estación terrestre para alertar al piloto al mando. En el caso de los drones que se encuentren bajo el modo de piloto automático, esta función es autónoma, y por lo tanto, deberá determinar y ejecutar por sí misma la maniobra apropiada para evitar colisiones. Además, también es necesario que sea capaz de retornar a la ruta original establecida, por lo que el nivel de control automático implementado en estos sistemas es significativo. 46 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Figura 3.12: Diferentes tipos de cámaras de visión artificial. U-Camera de Airelectronics, Hero 3 de GoPro, Tetracam ADC Light, Tau de FLIR (Grupo Acre) En la Figura 3.12 se pueden observar cuatro cámaras. Las dos primeras proporcionan imágenes dentro del espectro visible tienen peso y tamaños reducidos y gran calidad de imagen. La tercera se trata de una cámara multiespectral. Se utiliza sobre todo en cultivos, puesto que las imágenes infrarrojas (basadas en el principio explicado anteriormente) indican cambios en la vegetación antes de que aparezcan en el espectro visible. La última cámara se trata de una termográfica, basada en el principio de la radiación del cuerpo negro (en forma también infrarroja) en función de la temperatura. Los cuerpos con mayor temperatura emitan más radiación infrarroja que los que poseen menos temperatura. [12] 3.4. Sensores de luz Estos sensores son capaces de cuantificar la presencia de luz usando una serie de dispositivos como pueden ser las células fotoeléctricas. Son utilizados para medir la intensidad de luz incidente, y algunos de ellos son capaces de orientarse para mejorar el aprovechamiento de los rayos solares. Una aplicación práctica que se le podría dar, sería la siguiente: en los drones que utilizan la energía solar como fuente de energía, instalar este tipo de sensores podría optimizar la orientación de las placas solares para sacarle el mayor partido a la incidencia de los rayos del sol. Como se ha explicado anteriormente, dentro de este apartado se podrían incluir aquellas cámaras de vídeo con circuitería compleja, pero se decidió incluirla en el apartado de sensores de distancia y proximidad debido a su funcionalidad. Aquí se explicarán los sensores de luz sencillos, como pueden ser los fotodiodos, fototransistores y fotorresistencias. 1. Fotorresistencias. También llamado fotorresistor o LDR (del inglés, Light-Dependent Resistor), consiste en un dispositivo con una resistencia eléctrica cuyo valor varía en función de la luz incidente sobre él. Así, en los momentos en los que mucha (poca) luz incida sobre el sensor, esta resistencia disminuirá (aumentará). Esta variación suele producirse de manera no lineal, es relativamente lenta y además no se comporta de la misma forma al pasar de oscuro a claro que de claro a oscuro. Ofrecen mayor sensibilidad a la luz que los fototransistores. Típicamente se fabrican con un cristal semiconductor fotosensible, como podría ser el sulfuro de cadmio (CdS), puesto que son sensibles a un amplio rango del espectro visible y no visible (infrarrojos y ultravioleta). Figura 3.13: Esquema fotorresistencia 2. Fotodiodos. Este dispositivo de luz es comúnmente utilizado debido a que su respuesta a la luz es bastante rápida, de forma que se podría obtener un rango de tensiones lineal para un rango de 47 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados luminosidad. El fotodiodo está formado por un diodo semiconductor, construido como el diodo de unión-PN. En este caso, el semiconductor se encuentra expuesto a la luz debido a una lente transparente, y tendrá sensibilidad suficiente para la luz visible e infrarroja. El fundamento de un diodo se basa en tener un sentido normal de circulación de la corriente, llamado polarización directa. En el sentido contrario no permite pasar corriente. Pero en este caso, el fotodiodo se encuentra polarizado de forma inversa, es decir, la corriente que se modifica con los cambios de intensidad lumínica, circula en sentido inverso al permitido por la juntura del diodo. Esta aumenta cuando el sensor es excitado por la luz. Dependiendo de la aplicación se utilizarán unos fotodiodos u otros, cuyas diferencias residen en el material semiconductor utilizado (silicio, germanio, indio-galio-arsénico, sulfuro de plomo, etc.), con sus correspondientes rango de espectro. Figura 3.14: Esquema fotodiodo 3. Fototransistores. Al igual que el funcionamiento del fotodiodo se parecía al del diodo, este sensor funciona de manera similar a un transistor, con sus tres conexiones externas base-colectoremisor. De nuevo en este dispositivo se genera una corriente (colector-emisor) proporcional a la luz incidente en él (base-colector), gracias a una cápsula con una ventana transparente para el paso de la luz. De hecho, un transistor se puede convertir fácilmente en un sensor fototransistor conectando un fotodiodo entre colector y base. Pueden proporcionar una corriente mucho mayor que las de un fotodiodo estándar, ya que son más sensibles a la luz que los anteriores. Tienen, al igual que los fotodiodos, un tiempo de respuesta muy corto. Una característica importante de estos fototransistores es que proporciona variaciones mayores de corriente como respuesta a variaciones de intensidad luminosa, debido a que cuentan con un factor de amplificación. Figura 3.15: Esquema fototransistor Aplicaciones inmediatas que un UAV puede extraer de este tipo de sensores sería por ejemplo en almacenamientos de naves industriales. El vehículo con un sensor de luz, podría detectar cuando una serie de objetos apilados alcanzan una altura determinada para dar la orden de comenzar a amontonar en otra columna diferente. También se podría utilizar para conocer y rellenar huecos vacíos entre una serie de objetos en una estantería, lo cual ahorraría el tiempo que emplearía un ser humano en detectar dicho espacio vacío. 3.5. Sensores de Velocidad En el caso de los vehículos aéreos no tripulados, es necesario tener un control preciso sobre la velocidad a la que se está desplazando, para así poder mantenerla constante, anularla para dejarlo 48 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados estático o llevarla al máximo, por ejemplo. Se necesita que el sistema que acciona la velocidad de giro de los motores de sus hélices sea controlable desde un circuito remoto. En este apartado se hablará de los tacómetros otacogeneradores, los cuales son unos sensores capaces de medir velocidades angulares. Un tacogenerador sencillo se describe en la Figura 3.16, donde se observa que se utiliza un interruptor que utiliza fuerzas magnéticas para activarse o no. Esta activación dependerá de la cercanía o lejanía a la que se encuentre un diente magnetizado de la rueda dentada, de forma que cada vez que se acerque lo suficiente será accionado, y de esta forma se podrá obtener la velocidad de giro. Por cada vuelta, el dispositivo se activa y generará un voltaje de salida, el cual se utilizará posteriormente para calcular dicha velocidad. Figura 3.16: Tacogenerador: rueda dentada acciona el interruptor magnético [22] El tipo más utilizado en la actualidad de sensor de velocidad es el tacómetro de corriente continua. Funcionan de manera inversa a los motores de corriente continua: a partir de una velocidad angular son capaces de obtener una tensión de salida proporcionales a la velocidad de giro, es decir, consigue energía eléctrica a partir de mecánica. La Figura 3.17 explica el funcionamiento de este dispositivo. El rotor suele tener unas nbobinas, por lo que se producen 2ncontactos al haber 2 estátor. Si el eje gira con una velocidad de ω, se generará una tensión alterna con un valor de Vs=BfSsin (ωt)(3.5) donde Bfes la inducción del campo magnético y Sla superficie de la bobina en contacto. Después de calcular esta tensión alterna (proporcional a la velocidad de giro), se rectifica la tensión de manera que se obtenga una onda prácticamente constante, Vsen la Figura 3.17 Figura 3.17: Esquema del funcionamiento de un tacómetro. [22] La empresa SkyRC ofrece un modelo de tacómetro bastante apropiado para el aeromodelismo tanto de drones como helicópteros, el cual se puede observar en la Figura 3.18. 49 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Figura 3.27: Medida de los pseudorangos de 4 satélites usados para triangular la posición del receptor [3] En la Figura 3.27 se observan los 4 satélites necesarios para triangular la posición en 3D del receptor, en este caso un UAV. ¿Por qué son necesarias cuatro medidas del pseudorango? Pues bien, para calcular una posición en tres dimensiones bastaría con utilizar tres medidas diferentes, así se obtendrían las coordenadas (x,y,z) o(latitud, longitud, altitud). Pero para resolver el problema del offset de los relojes, se necesita una medida adicional. De esta forma se obtendría un sistema de cuatro ecuaciones no lineales con cuatro incógnitas: latitud, longitud, altitud y el offset del reloj del receptor. El GPS es el más desarrollado y usado de los diferentes sistemas GNSS (del inglés, Global Navigation Satellite System). En Febrero del 2012, el Congreso de los Estados Unidos instó a la Agencia Federal de Aviación FAA, (del inglés Federal Aviation Agency), que llevara a cabo un plan para la seguridad de la rápida incorporación de los UAVs dentro del ámbito civil. Se pedía un aumento de seguridad para la navegación de este tipo de vehículos puesto que utilizaban la señal de GPS de uso comercial y no la de uso militar, la cual tiene una mayor precisión y seguridad. Como ejemplo, la Universidad de Texas demostró que era capaz de acceder a la señal GPS de uso civil y modificarla, introduciendo información falsa en dichas señales (fueron capaces de desviar el curso de un velero cientos de metros que se encontraba en el Mediterráneo). Una medida para no estar completamente expuestos a estos ataques de intrusismo, tal y como se hace actualmente y se propone en [19], es utilizar fuentes independientes que proporcionen información sobre la navegación. Por ejemplo , el avión comercial B-787 actualmente utiliza varios sistemas de GNSS como GPS (EEUU), BeiDou (China), GLONASS (Rusia) y Galileo (Europa). Con varias señales que proporcionan la misma información es más complejo acceder y modificar los datos. Al igual que en todos los sistemas sensoriales, el GPS no se libra de poseer fuentes de errores que son necesarias tener en cuenta para mitigar sus efectos. La precisión de la medida que se obtiene de la posición dependerá tanto de la precisión de los pseudorangos como de la geometría de los satélites que el receptor tiene a su alcance. La mejor posición es formando un tetraedro con el receptor en el centro de una de las caras, y el vértice opuesto a dicha cara sobre el usuario. Para tener en cuenta los errores debidos a la geometría, se utilizan una serie de factores denominados DOP (del inglés, Dilution of Precision). También se producen errores debido al tiempo de propagación de la señal. Puesto que estas viajan a la velocidad de la luz, porque son señales electromagnéticas, un error de 10ns podría resultar en un error en posición de hasta 3m. Algunas fuentes que provocan errores en los tiempos de propagación son brevemente descritos a continuación. Datos de la efemérides. La efemérides es una descripción matemática de la órbita del satélite. Para calcular la posición del receptor es necesario conocer la del emisor con precisión, es decir, no cometer errores en esa descripción matemática sobre su órbita. Reloj del satélite. Los relojes utilizados son atómicos, fabricados de cesio y rubidio. A lo largo del día arrastran un error de 10ns, por lo que como se dijo anteriormente provocan alrededor de 3m de error en posición. Este problema se soluciona actualizando los relojes cada 12 horas, por lo que el error en posición medio es de 1 o 2 metros. 56 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Ionosfera y Troposfera. La ionosfera es la capa más externa de la atmósfera, donde se producen los mayores errores en la propagación (alrededor de 2 y 5 metros) debido a la presencia de electrones libres. La troposfera, por el contrario, es la capa más baja de la atmósfera, donde se concentra la mayor masa de la atmósfera. Ahí se produce la mayor actividad meteorológica, por lo que los cambios de presión, temperatura y humedad afecta a la velocidad de la señal enviada. Los errores que ésta introduce son cercanos al metro. Multitrayectoria. Se produce cuando el receptor recibe señales reflejadas en edificios u otros grandes obstáculos que provocan errores por debajo de un metro en la mayoría de las circunstancias. Estas señales reflejadas ocultan y dañan la medida real de la posición. La cantidad de información que se obtiene a partir del uso de GPS es bastante llamativa. En el propio código de la señal que llega al receptor se indican valores no solo de posición, sino también de velocidad, número de satélites utilizados, hora UTC (del inglés, Coordinated Universal Time), DOPs, elevación y azimut de satélites, calidad del mensaje, etc. Todos estos apartados se verán con más profundidad a la hora de analizar los datos obtenidos en el formato NMEA proporcionado por el receptor GPS. 3.9.1. Segmentos Este sistema de posicionamiento y navegación consiste en tres segmentos: segmento espacial, segmento de control y segmento de usuario. Las Fuerzas Aéreas de los Estados Unidos son las encargadas del desarrollo, mantenimiento y operación de las dos primeras. Segmento Espacial. Se trata de la constelación de los satélites transmitiendo radio señales a los usuarios. Consta de 24 satélites operacionales el 95% del tiempo. Se encuentran distribuidos en 6 planos orbitales, con 4 satélites situados en cada plano. Esta disposición permite que cada usuario pueda ver al menos cuatro satélites desde cualquier punto del planeta. Vuelan en órbitas medias denominadas MEO (del inglés, Medium Earth Orbit), las cuales son circulares, y a una altura de 20.200km aproximadamente. Cada satélite orbita alrededor de la Tierra dos veces al día. Segmento de Control. Consiste en la red de instalaciones que monitorizan los satélites de GPS, interpretan las transmisiones, realizan análisis y mandan instrucciones y datos a las constelaciones. Son las encargadas de actualizar la posición real de los satélites (las efemérides) basándose en observaciones, además de sincronizar los relojes atómicos. Las estaciones de control se encuentran distribuidas por todo el planeta, y consta de: una estación de control principal (en Colorado, EEUU), una estación de control principal alternativa (en California, EEUU), 11 antenas de control y orden, y 15 lugares de monitorización. 57 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Vandenberg AFB California Alternate Master Control Station Air Force Monitor Station Hawaii Master Control Station Schriever AFB Colorado NGA Monitor Station South Korea Australia Bahrain South Africa United Kingdom Ecuador USNO Washington Alaska New Zealand AFSCN Remote Tracking Station Ascension Diego Garcia Cape Canaveral Florida Kwajalein Ground Antenna New Hampshire Greenland Guam Updated April 2016 GPS Control Segment Figura 3.28: Segmentos de Control de GPS Segmento de Usuario. El GPS se ha convertido en los últimos años en uno de los pilares tecnológicos en los que se basan numerosas aplicaciones que se utilizan de forma cotidiana. El segmento de usuario consta de los receptores que poseen los usuarios de GPS para obtener su posición a partir de las señales recibidas. Este receptor debe contener un reloj de cuarzo para sincronizarlo con los relojes de los satélites y así obtener una medida más fiable de la posición. Estos receptores deben también decodificar el mensaje enviado por el satélite. Suelen estar formados por: antenas, receptores, microprocesadores, almacenamiento de datos para el cálculo de rutas, unidades de control para facilitar el uso al usuario, entre otras. Algunas de las prestaciones de este dispositivo son las siguientes: •Actualización de la posición de 0.5 a 1 segundo •Precisión en torno a los 15m •Precisión de 0.1m/s en cuanto a la velocidad del usuario •Precisión de 100ns en cuanto a la referencia temporal 3.9.2. Disponibilidad, Integridad y Continuidad La principal limitación con la que cuenta este sistema es con la dependencia de los Estados Unidos, más en concreto con el departamento de defensa. Se están desarrollando otros sistemas de posicionamiento global como son el Galileo (Europa) o GLONASS (Rusia) para paliar esta dependencia, ya que en cualquier momento podrían decidir limitar el uso del GPS por parte de la población civil. Actualmente el sistema GPS es el que consigue una mayor precisión con respecto al resto de sistemas de navegación en uso. Sin embargo, la precisión de un sistema no es la única variable a tener en cuenta a la hora de utilizar un sistema u otro, puesto que existen otros conceptos como pueden ser su disponibilidad, integridad y continuidad. En este apartado se describen brevemente. Disponibilidad. Porcentaje de tiempo que dicho sistema se encuentra utilizable, dentro de su área especificada de cobertura. El concepto de utilizable hace referencia a que cumpla unos requisitos mínimos previamente especificados para que la información que proporciona sea mínimamente fiable. Un ejemplo típico sería limitar alguno de los DOP a valor concreto, por debajo del cual no se tendrá en cuenta la señal recibida y por tanto los datos contenidos en ella. El área de cobertura del GPS incluye toda la superficie terrestre, aunque hay que tener en cuenta el ángulo de elevación del satélite en el horizonte a partir del cual se considera visible para el receptor, llamado “ángulo de máscara”. Éste será mayor en entornos urbanos que despejados. 58 3 Modelado de Sensores en los Vehículos Aéreos no Tripulados Integridad. Capacidad del sistema para controlar que el propio sistema no debe ser utilizado, ya sea porque no se encuentre operativo o porque posea errores. Se le considera un parámetro que equivaldría a una medida de confianza que se pueda tener en el sistema. Los efectos de radiación del espacio que alteren los relojes o la electrónica, los fallos de los satélites, o el error humano de software o hardware en el segmento de control son algunas de las causas que provocan que la integridad del sistema no se la esperada. La probabilidad de que sucedan este tipo de fallos son bajos, pero inadmisibles para aplicaciones de navegación aérea. Puesto que se trata de un parámetro realmente importante puesto que en él reside la confianza en el sistema, se han desarrollado diversos procedimientos para aumentar la integridad del GPS. Una de las técnicas más conocidas es el GPS diferencial oDGPS (del inglés, Differential GPS). Se trata de utilizar una serie de estaciones en tierra (GBAS, del inglés, Ground-Based Augmentation Systems) cuya posición se conoce con gran precisión. Estas están equipadas con un receptor GPS y se encuentran en comunicación con el usuario, para proporcionar mayor precisión a la señal del satélite. Figura 3.29: Ejemplo de DGPS basado en GBAS Continuidad. Probabilidad de que el sistema de navegación pueda ser usado de manera continua durante la realización de una misión u operación. Que pueda ser usado hace de nuevo referencia al concepto previamente definido de “utilizable”, es decir, que cumpla unos requisitos mínimos impuesto por la operación a realizar, y que de nuevo podrían ser expresados en términos de los DOPs. Dependiendo de la misión u operación que se lleve a cabo, la continuidad será una u otra, pero en cualquier caso se encontrará íntimamente ligada con fallos no planificados en los satélites. La probabilidad de que en un satélite se produzca un fallo que produzca que dicho satélite deje de emitir de forma no planificada es del 0.0001 %. Muchos sistemas militares de navegación inercial podrían ser reemplazados con sistemas con menos fiabilidad y coste si se garantizara la continua disponibilidad del GPS. Esta continua disponibilidad serviría para actualizar al sistema inercial y limitar así su error, el cual se propaga a lo largo del tiempo. Puesto que no está garantizada, se desarrollan sistemas alternativos de navegación inercial más costosos pero cuya fiabilidad es mucho mayor, debido a que se conoce que no va a sufrir ataques externos. [19] 59 4|Filtro de Kalman Esta sección se centrará en la explicación, formulación y análisis de una herramienta muy comúnmente utilizada, el Filtro de Kalman. Será necesario comprender sus expresiones y algoritmos para el problema que se encuentra bajo estudio. Es muy típico encontrar este tipo de procedimiento en multitud de sistemas, ya que ofrece una amplia gama de aplicaciones en diversos ámbitos. En este caso, puesto que el objetivo final es implementar en Matlab una simulación mediante la cual se puedan integrar dos señales recibidas por un UAV, señal GPS y señal de una IMU. Se explicarán también otros tipos de filtros utilizados actualmente. (EXTENDER CUANDO TERMINES) 4.1. Introducción El filtro de Kalman fue desarrollado en primer lugar para su uso en la navegación de aeronaves, pero su implementación se ha visto ampliada, abarcando todo tipo de ámbitos y campos. Su principal uso reside en estimar los estados de un sistema, de los cuales se ha hablado con anterioridad en este trabajo. Dichos estados suelen tener una propiedad común la cual hace necesaria que se le apliquen este tipo de filtros: solamente pueden ser observados de forma imprecisa por el propio sistema. Esto es debido a que el estado de un sistema es proporcionado por los datos procedentes de los sistemas sensoriales anteriormente detallados, los cuales se han visto que poseen diversas fuentes de ruido que producen errores en las medidas. Por este motivo se ha decidido aplicar este tipo de sistema de filtrado para mejorar la calidad de las señales recibidas (en el caso de GPS) o de los datos obtenidos (en el caso de la IMU). Existen otro tipo de razones por las cuales ha sido elegido este filtro, como su atractivo debido a que minimiza la varianza del error de la estimación, las cuales serán explicadas más adelante. El objetivo de este apartado es introducir los conceptos básicos sobre los algoritmos del Filtro de Kalman, y utilizarlos en una simulación para resolver un problema real sobre la navegación de un vehículo, en este caso aéreo y no tripulado. Este método proporciona una herramienta suficientemente fiable para estimar la posición (o cualquier otra variable a controlar) actual del vehículo. En el momento de la obtención de datos mediante las señales mencionadas anteriormente, es inevitable que se produzcan efectos indeseados de ruido, los cuales pueden corromper el mensaje. Un buen algoritmo de filtrado debería ser capaz de eliminar este ruido de las señales (el acumulado al utilizar un acelerómetro, al atravesar la atmósfera la señal de GPS, etc.) a la vez que retiene la información útil que se desea obtener. Es por ello que anteriormente se han detallado las fuentes de errores de todo tipo de sensores, porque un mejor conocimiento del origen de este error proporciona información necesaria para diseñar un modelo lo más realista posible. 4.2. ¿Qué es el filtro de Kalman? El filtro de Kalman es un proceso matemático iterativo que utiliza una serie de ecuaciones y datos de entrada constantes para rápidamente estimar el valor real de la posición, velocidad, o cualquiera que sea la variable, del objeto medido, cuando los valores medidos contienen errores impredecibles o aleatorios, incertidumbre o variación. Se dice que es un proceso rápido puesto que es capaz de acercarse al valor real tomando un número relativamente bajo de esos datos de entrada. 60 4 Filtro de Kalman Figura 4.1: Proceso iterativo del Filtro de Kalman En la Figura 4.1, obtenida de [24], se muestra a grandes rasgos, un diagrama de los principios en los que se basa en filtro de Kalman. Será realmente útil cuando llegue el momento de explicar las tres principales ecuaciones o cálculos necesarios. En la misma figura se puede ver que de hecho se trata de un proceso iterativo. El primer término que se necesita calcular, en cada iteración, es el llamado “Ganancia de Kalman”. Después, será calculado la “Estimación Actual”, la cual significa que se procederá a actualizar el valor estimado. Por último, el “Nuevo Error en la Estimación” será recalculado a partir de los dos cálculos anteriores. Se utilizará por comodidad el término error siempre como sinónimo de incertidumbre. Los cálculos de estos tres términos determinarán las tres principales ecuaciones mencionadas anteriormente. Estos son los requisitos y la funcionalidad de cada uno de los términos: Ganancia de Kalman. Se necesitarán dos datos, el error en la estimación y el error en la medida, como se puede observar en el diagrama de la Figura 4.1 (dos flechas llegan a la ganancia de Kalman). •Primero, el error en la estimación, que puede ser tanto el error en el anterior o si es el caso de la primera iteración, el error original. La estimación, como se verá más adelante, está basada en el modelo dinámico que se emplee para modelar el sistema. Por ejemplo, en el caso del seguimiento de una trayectoria de un UAV, el modelo se regirá por las leyes de Newton de la Cinemática x=x0+vt + 1/2at2. El error en la estimación se produce en el intento de pasar a modelo matemático un problema real, puesto que se cometen fallos tanto debido a aproximaciones, como a redondeos, simplificaciones, etc. •Segundo, el error en el dato de entrada recibido en ese instante. Éste es debido a los errores que se cometen al tomar medidas con sensores, y dependerá del sistema sensorial que se esté utilizando. En el apartado anterior se han explicado los diferentes errores que se cometen con los diferentes sensores. En el caso bajo estudio, donde se tomarán datos de un GPS y una IMU, los errores serán los correspondientes a dichos sistemas. Por comentar algún ejemplo, el GPS comete errores de sincronización entre sus relojes, y la IMU acumula errores a lo largo de la trayectoria. El objetivo de la Ganancia de Kalman será ponderar una importancia relativa en el primer error (de la estimación EEST ) versus el segundo error (dato de entrada), llamado en este caso EMEA 61 4 Filtro de Kalman debido al inglés measurement. KG =EEST EEST +EMEA KG ∈[0,1] (4.1) Se ha querido explicar el significado de esta ecuación mediante el esquema de la Figura 4.2. La Ganancia de Kalman puede tomar valores comprendidos entre 0 (baja) y 1 (alta). ¿Qué significa obtener una ganancia de Kalman alta? Como se observa en el esquema, una ganancia de Kalman alta implica que el Error de la Estimación es grande comparado con el Error de la Medida. Y si el Error de la Medida es más pequeño que el Error que se ha estimado, interesa tener más en cuenta la medida (dato obtenido), y “despreciar” el valor de la estimación porque contiene un error alto. Este concepto de dar más importancia a un valor que a otro es la función principal de la ganancia de Kalman, como se verá al realizar el siguiente cálculo: Estimación Actual. Figura 4.2: Explicación del uso de la Ganancia de Kalman Estimación Actual. Para el cálculo de la estimación actual ESTtse necesitarán tres variables (tres flechas llegan a la estimación actual), tal y como se puede ver en el diagrama de la Figura 4.1: la estimación anterior ESTt−1, el valor medido por el sensor MEA, y la ganancia de Kalman calculada anteriormente KG. Esta estimación se rige por la siguiente expresión, en las cuales aparecen dichas variables ESTt=ESTt−1+KG[MEA −ESTt−1] =ESTt−1[1 −KG] + KG ∗MEA (4.2) •Primero, la estimación anterior. Como se ha dicho con anterioridad, el filtro de Kalman es un proceso iterativo, y por lo tanto, en cada iteración, cualquiera que sea la estimación actual para un tiempo tpasará a ser la estimación anterior en el tiempo t+1. En el diagrama de la Figura 4.1 se observa que la estimación anterior se puede obtener de dos formas: si es la primera iteración se le asignará un valor original que se considere oportuno, y si se trata de la segunda o siguientes iteraciones, corresponderá a la anterior estimación actual. Una de las características más importantes del filtro de Kalman que lo hacen tan especial reside en dicho “valor original” que se le asigna al comienzo del proceso. Este filtro tiene la propiedad de que no importa qué valor es determinado como inicial, puesto que rápidamente será capaz de acercarse al valor real basándose en los errores. •Segundo, el valor medido por el sensor. Este es el dato de entrada que se recibe de los sensores, en este caso, del GPS y/o de la IMU. Estos sensores proporcionarán datos sobre 62 4 Filtro de Kalman posición y velocidad que, como ya se ha explicado, contienen fuentes de errores, y según el grado de error que contengan las señales, se tendrán más o menos en cuenta a la hora de calcular la nueva estimación. •Y tercero, la ganancia de Kalman. Pues bien, ha llegado el momento de decidir. En este punto se tiene: una estimación anterior (a partir de la cual, con el modelo matemático correspondiente se conseguirá calcular la estimación actual) con su correspondiente error, y unos nuevos datos de entrada con sus correspondientes errores. ¿Cómo decidir a qué valor de estos darle más importancia? ¿Se toma uno y se desprecia el otro por completo? Para responder a estas preguntas se dispone de la ganancia de Kalman y de la ecuación 4.2. La ganancia de Kalman será la herramienta utilizada para ponderar entre ambos valores basándose en cuál de los dos posee un mayor error, como se indicó en la Figura 4.2. Atendiendo a la ecuación, se observa cómo, ◦si la ganancia de Kalman es alta (KG=1, error en la estimación alto), el término ESTt−1[1 −KG]se anularía puesto que contiene alto error, mientras que el término KG ∗MEA se mantendría porque el error en la medida es bajo (en comparación con el de la estimación). ◦si la ganancia de Kalman es baja (KG=0, error en la medida alto), el término KG ∗ MEA se anularía porque contiene alto error, mientras que el término ESTt−1[1 −KG] se mantendría puesto que el error en la estimación es bajo (en comparación con el de la medida). Nuevo Error en la Estimación. Para este apartado se necesitarán dos variables (dos flechas llegan a este punto en la Figura 4.1): la estimación actual y la ganancia de Kalman. Es decir, los dos pasos anteriores son necesarios para calcular este error. EESTt=EMEAEESTt−1 EMEA +EESTt−1 =⇒EESTt= [1 −KG]EESTt−1(4.3) Es en este punto en el que comienza la nueva iteración, y por lo tanto el término EESTt−1, hace referencia al error en la estimación calculada justo en el apartado anterior, en donde se hablaba de actual pero que ya ha pasado a ser la anterior. Es interesante destacar que el error en la estimación actual siempre será más pequeño que el anterior, debido a que KG ∈[0,1]. Cuanto mayor sea KG, más rápido encontrará el valor verdadero. En cada iteración, algunos valores son obtenidos del proceso iterativo, los cuales serán usados para obtener un mejor resultado de aquello que se esté calculando. En el caso bajo estudio será la posición y la velocidad de un sistema UAV mediante simulación de la toma de datos de un GPS y una IMU. Cuantas más iteraciones se realicen, más precisos serán estos valores obtenidos. 4.3. El modelo multidimensional En la sección anterior se introducen unas breves pinceladas de sobre qué trata el filtro de Kalman, así como las ecuaciones y cálculos básicos con el fin de entender el proceso. Pero el problema bajo estudio es más complicado que eso, puesto que es necesario tener en cuenta no sólo una medida en cada instante de tiempo, sino que se obtienen varios datos de entrada (tres medidas para las tres direcciones de la posición, otras tres para las de la velocidad, etc.). Por ello, es necesario explicar el modelo multidimensional de este tipo de filtros. 63 4 Filtro de Kalman Figura 4.3: Proceso Iterativo que explica el modelo multidimensional mediante las ecuaciones del algoritmo En la Figura 4.3 se intenta representar los diferentes pasos a seguir a la hora de implementar el filtro de Kalman mediante las ecuaciones que componen el algoritmo. Se intentará desarrollar cada uno de los recuadros señalados indicando qué significa cada variable y para posteriormente aclarar cómo se calcula cada una de ellas. a) Estado Inicial. Estado desde el cual parte la simulación o el modelado del sistema. Consta de x(0) = Matriz de estados iniciales. En dicha matriz se incluirán los valores iniciales de las variables que se quieren filtrar. En el caso del movimiento de un vehículo, normalmente se estudiarán la posición y la velocidad, y por lo tanto, esta matriz constará de 6 elementos: las 3 coordenadas en posición y las 3 coordenadas en velocidad. P(0) = Matriz de la covarianza inicial del proceso. Hace referencia a los errores en las estimaciones, de los cuales se habló anteriormente. Este es el término que corresponde con el numerador de la ecuación 4.1. Posteriormente se explicará cómo se construye esta matriz. b) Estado Anterior. En la primera iteración el estado inicial se convierte en el estado anterior. Una vez que se comience a realizar el proceso completo, este estado anterior vendrá definido por los cálculos propios del estado en un instante (actual) y no del inicial. En este punto las matrices x(k-1) yP(k-1) significan lo mismo que en el paso anterior, matriz de estados y matriz de covarianza del proceso (error en la estimación). c) Nuevo Estado. A estas alturas ya se puede predecir el estado actual kpen el que se encuentra el sistema. Las ecuaciones correspondientes a este punto son las siguientes x(0) =⇒x(k−1) =⇒x(kp) = Ax(k−1) + Bu(k) + w(k) P(0) =⇒P(k−1) =⇒P(kp) = AP(k−1)AT+Q(k)(4.4) donde 64 4 Filtro de Kalman x(k) = Matriz de estados del sistema en el instante actual k(posición actual, velocidad actual, o cualquier otra variable que se esté considerando). Esta primera ecuación de la expresión 4.4 se denomina “Ecuación de estado”. u(k) = Matriz de variables de control en el instante actual k. En el caso bajo estudio, donde se pretende controlar un vehículo que vuela al aire libre, una de esas variables de control podría ser por ejemplo la gravedad. Esta aceleración de la gravedad es constante y afecta directamente a la posición del vehículo en cada instante, por lo que sería necesario modelarla como una variable de control. w(k) = Matriz del ruido en el estado calculado. Como se ha hablado anteriormente, pasar de un modelo real a un modelo matemático implica errores de aproximaciones, redondeos o simplificaciones que deben ser tenidos en cuenta en algún punto. De ello se encarga esta variable w(k). Cuanto más próximos estén el modelo simulado y el modelo real, menor será esta variable. P(k) = Matriz de la covarianza del proceso. La función es la misma que la explicada en el estado inicial: el error en la estimación calculada. Se verá cómo es necesario conocer la matriz Ay su traspuesta para actualizarla en cada iteración. Esta depende de la covarianza calculada anteriormente y de la covarianza del ruido del proceso Q A,B = Matrices utilizadas para convertir el estado de entrada en el estado nuevo (A), y las variables de entrada en estado nuevo (B). Dependiendo del modelo que se realice del sistema será necesario utilizar unas matrices u otras. Más adelante se explicará también cómo se forman estas matrices dependiendo de las variables que se toman como estados y del modelo utilizado. Q= Matriz de la covarianza del ruido del proceso. La función principal de Qes prevenir que la matriz Pde covarianza del estado se convierta en un valor muy pequeño o eventualmente que llegue a cero. Si Pllegase a valer cero, significaría que las medidas tomadas por los sensores que se están utilizando serían ignoradas por completo, aunque tuvieran una información de gran valor. [24] En la primera ecuación de 4.4, se puede observar, por tanto, cómo para calcular el nuevo estado (por ejemplo, la nueva posición del objeto), es necesario conocer el estado anterior (la posición anterior), la aceleración que sufre el vehículo (o cualquier otra variable de control que afecte a su trayectoria), y el error producido por modelar un sistema real (si se utilizan las leyes de la Cinemática se desprecia el rozamiento, por ejemplo). d) Entrada de la medida. Aquí entra en juego el papel de los sensores. Estos son los que toman las medidas y en este paso es en el que esas medidas se procesan para incluirlas en el proceso de filtrado de la señal. Para procesar esa medida es necesario tener en cuenta la siguiente expresión y(k) = Hx(km) + z(k)(4.5) que es la denominada “Ecuación de salida”, y donde y(k)= Es el valor que se introduce en el proceso de filtrado, el cual corresponde a la medida tomada por el sensor x(kM), pero a la cual se le añade el correspondiente valor del error que comete el sensor z(k). El nombre puede dar lugar a confusión, puesto que comúnmente se denomina ecuación de salida (salida del sensor, puesto que saca datos), pero que aquí se le ha nombrado como entrada de la medida (puesto que es el dato que entra en el filtro) x(kM)= Medida (de ahí el subíndice kM) tomada directamente por el sensor. Esta variable típicamente es un vector de datos, como por ejemplo, las 3 datos para la posición y 3 datos para la velocidad. 65 5 Caso práctico - Simulación en Matlab del Filtro de Kalman x. Este hecho será importante para explicar el resultado obtenido al finalizar la simulación del filtro. w= Error en el proceso. Q=Pn−1 i=0 (wi−¯w)(wi−¯w)T n−1=⇒donde ¯w=Pn−1 i=0 (xt−1−xt−i−1) n wi= (xt−1−xt−i−1) (5.2) R= Matriz de covarianza del ruido de los sensores. z= Error en la medida. R=Pn−1 i=0 (zi−¯z)(zi−¯z)T n−1=⇒donde ¯z=Pn−1 i=0 (yt−i−Ht−ixt−i) n zi= (yt−i−Ht−ixt−i) (5.3) Los cálculos realizados en Matlab se pueden encontrar en el apéndice B. Según [11] y [24], a la matriz de covarianza del proceso P se le puede asignar un valor cualquiera puesto que su valor se va actualizando en cada iteración del filtro y tiene la peculiar característica de que converge rápido a su valor real. También en esta literatura se indica que se utiliza una matriz diagonal, de forma que se imponga la no existencia de correlación entre las variables de estado. En el caso bajo estudio, esto significa que el desplazamiento producido en una variable de estado xno afecta al desplazamiento en las otras dos variables de estado yyaltitud, y con estas dos el mismo caso. Con esto, se ha decidido asignar a P el valor siguiente P=  700 070 007  (5.4) 5.3. Resultados obtenidos de la simulación Una vez que se han definido todos los parámetros para implementar el filtro de Kalman, se comienza con el proceso iterativo, tal y como se mostraba en el diagrama de la Figura 4.3. De esta forma, tal y como se puede vislumbrar en el código ejecutado en el apéndice B, se pueden obtener unos datos sobre una trayectoria estimada, basada en los datos obtenidos por el GPS y por Google Maps. El filtro de Kalman asegura obtener un resultado en el cual la varianza del error es la mínima, es decir, la trayectoria estimada una vez llevado a cabo el algoritmo del filtro será la más parecida a la trayectoria real (la que cuenta con menos error). La Figura 5.5 se ha girado y colocado de forma que se vea lo más claro posible el efecto del filtro sobre las señales. La gráfica en color verde discontinua se corresponde con los datos obtenidos del GPS. Como se puede observar, la altitud proporcionada por el GPS no es muy fiable, puesto que se conoce que la trayectoria calculada fue tomada en un puente, el cual no varía de altitud del orden de los 8 metros (de 166m a 174m) tal y como se ve en la figura. 72 5 Caso práctico - Simulación en Matlab del Filtro de Kalman La gráfica en color azul discontinua corresponde con los datos calculados a partir de Google Maps. Se observa claramente que se trata de un movimiento rectilineo, puesto que la altura es completamente constante e igual a 167m. Ésta difiere un par de metros de los 169m que se obtienen de media si se toman los datos del GPS. Esta trayectoria se puede decir que es aquella que se intentaba obtener cuando se tomaron las muestras. Se pretendía circular a velocidad constante en un tramo donde la altura con respecto al nivel del mar fuese lo más constante posible, por lo que surgió la idea de tomar muestras sobre un puente. ¿Cuál es el problema? Que existen múltiples factores que introducen error en estos datos. Para empezar, es muy difícil controlar que la velocidad de un vehículo sea completamente constante. Por ello se decidió utilizar un coche que tuviese la opción de conducir con velocidad de crucero, pero aun así no se asegura que dicha velocidad no varíe en ningún punto. A la hora de controlar el coche, también una fuente de error reside en los pequeños desplazamientos que se producen por movimientos en el volante. También es posible que se desvíe de la trayectoria ideal debido a fuerzas externas, como pueden ser la acción del viento o irregularidades en el terreno. Y por supuesto, los valores han sido tomados por una herramienta como Google Maps, la cual no se libra de errores a la hora de proporcionar la posición del punto que se señale en el mapa. Con esto se pretende decir que esta “señal ficticia” dista también, aunque en menor medida, de la trayectoria real seguida por el vehículo. 4.93454.9354.9355 x 106 −6.046−6.0455−6.045−6.0445−6.044−6.0435−6.043−6.0425 x 105 165 166 167 168 169 170 171 172 173 174 Eje Y Filtrado de señales de GPS y Google Maps Eje Z Datos GPS Google Maps Kalman Figura 5.5: Comparación de las señales GPS, Google Maps y filtro de Kalman en altitud Obtener una trayectoria que diste lo menos posible de la trayectoria real es el objetivo por el que se ha realizado esta simulación. Intuitivamente, la trayectoria real que se ha seguido con el vehículo será bastante parecida a la trayectoria calculada por Google Maps, en el sentido de que sería muy similar a una línea recta, puesto que las desviaciones producidas por los factores de error explicados en el apartado anterior son pequeños en este caso.En el caso de estudiar un UAV de verdad, dichos factores cobrarían más importancia debido a que son más sensibles a perturbaciones externas, debido a su menor peso, tamaño y estabilidad. Es decir, es más fácil que se desvíe un UAV 3 metros debido a ráfagas de aire que a un coche. Dicha trayectoria casi rectilínea sería interesante dibujarla también en la Figura 5.5, para de verdad comprobar lo que se pretende en esta simulación: que la trayectoria de Kalman es la más acertada a la hora de estimar una trayectoria real. Centrados de nuevo en la Figura 5.5, se observa cómo la trayectoria estimada por el filtro de Kalman (en color rojo, continua) se encuentra en todo momento en un punto intermedio entre la señal de GPS y la señal de Google Maps. Se recuerda que esto es debido a que la ganancia de Kalman era un factor que ponderaba la importancia relativa que se le debe dar a una señal con respecto a la otra, basándose en el error en cada momento. Es decir, en cada iteración cuando se obtienen nuevos 73 5 Caso práctico - Simulación en Matlab del Filtro de Kalman datos, la ganancia de Kalman hará que el nuevo dato de la estimación conste de un porcentaje de una señal (por ejemplo, ganancia de Kalman=0.6 implica que tomará un 60 % de la señal GPS) y de otro tanto por ciento de la otra (en este caso tendría un 40% de la señal de Google Maps). Se observa cómo al principio de la gráfica la trayectoria de Kalman se encuentra en un valor más cercano al GPS y luego se aproxima cada vez más a la señal de Google Maps. El inicio de la trayectoria de Kalman es tan inexacto debido a que se definió en primer lugar una P aleatoria y de valor alto para, de hecho, comprobar que dicha matriz converge y se va acercando cada vez más al valor que posee menor error, en este caso la trayectoria de Google Maps. Para visualizar mejor este hecho, la Figura 5.6 se corresponde con una matriz P que en lugar de ser una diagonal de valor 7, es una diagonal de valor 1. Se observa cómo la gráfica de la trayectoria de Kalman comienza más próxima al valor más acertado, en este caso la señal de Google Maps, la cual posee menos errores que la de GPS (en altitud). 4.93464.93474.93484.93494.9354.9351 x 106 −6.046−6.0455−6.045−6.0445−6.044−6.0435−6.043−6.0425 x 105 165 166 167 168 169 170 171 172 173 174 Eje Y Filtrado de señales de GPS y Google Maps Eje X Eje Z Datos GPS Google Maps Kalman Figura 5.6: Comparación de las señales GPS, Google Maps y filtro de Kalman en altitud, P más baja Anteriormente se dijo que sería interesante el valor de la matriz Q para interpretar el resultado, y aquí se ve claramente por qué: se observa cómo rápidamente el filtro decide que debe acercarse más a la trayectoria de Google Maps, por lo tanto, considera que dicha medida posee menor error que la medida del GPS. Puesto que los errores residen en las matrices Q y R, habrá que analizarlas para encontrar el por qué de este fenómeno. Pues bien, comparando ambas, se observa que la matriz Q tiene la siguiente forma Q=  **0 **0 000  (5.5) Este fenómeno se debe a que los errores dependen de cuánto varían las medidas con respecto a su valor medio (así se calcula la covarianza, ecuación 5.2), y se sabe que la altitud calculada según Google Maps es constante, por lo que dicha diferencia será nula. Este hecho de tener un valor nulo hace que la ganancia de Kalman siempre se vaya a decantar un poco más por este valor que por otro que sí tenga dicho error, por muy pequeño que sea será mayor que éste. Parece que, puesto que el GPS proporciona medidas bastante erróneas de altitud, que no se debería confiar en el valor proporcionado por dicho sistema. Pero en la Figura 5.7, se ha intentao orientar la figura para que se vea en los ejes (x,y), tal y como se veía en la imagen de Google Earth 5.3. 74 5 Caso práctico - Simulación en Matlab del Filtro de Kalman 4.9347 4.9347 4.9348 4.9348 4.9349 4.9349 4.9349 4.935 4.9351 4.9351 x 106 −6.046 −6.0455 −6.045 −6.0445 −6.044 −6.0435 −6.043 −6.0425 x 105 Filtrado de señales de GPS y Google Maps Eje Y Eje X Datos GPS Google Maps Kalman Figura 5.7: Comparación de las señales GPS, Google Maps y filtro de Kalman en ejes (x,y) En esta figura ya es más difícil diferenciar qué señal se encuentra más próxima a la trayectoria estimada por el filtro de Kalman. También parece ser que desde esta perspectiva se aprecia que sigue mejor a la señal de Google Maps, pero en este caso no dista tanto de la señal de GPS. En este caso sí que se podría afirmar que GPS proporciona datos bastante fiables, acertados con la posición “real”, en el caso de que la señal de Google Maps fuera completamente la real. En cualquier caso, si se quisiera conocer con mayor exactitud el camino real seguido por el vehículo, habría que mirar la gráfica de color rojo, la cual se corresponde con la trayectoria estimada por el filtro de Kalman, puesto que es la que menos errores posee. 75 6|Conclusiones y Líneas Futuras Con este estudio se ha intentado analizar el estado actual en el que se encuentran los sensores y los sistemas que integran un UAV. Para comenzar el proyecto, se examinó la herramienta de espacio estados, puesto que se considera el instrumento más útil a la hora de trabajar con sistemas multivariables, como pueden ser los vehículos aéreos no tripulados. Con el estudio sobre el espacio de estados, se obtiene una visión general de cómo plantear un problema real de seguimiento de trayectorias, es decir, qué variables se deben tomar como estados puesto que definen el sistema modelado. Para ello, fue necesario introducir un apartado sobre los diferentes tipos de sensores que se encuentran implementados en los drones. Se obtiene de este apartado que el desarrollo de esta tecnología tan puntera, no depende sólo de un ámbito de la ingeniería como puede ser la aeroespacial, sino que se trata de un entorno en el que se necesita trabajar con grupos de ingenieros de muchas especialidades, como pueden ser automáticos, electrónicos, de aviónica, mecánicos, de energía, etc. Al estudiar dichos sistemas, se observó que este sector de los UAVs está bastante avanzado en el ámbito militar, como es lógico, puesto que allí se desarrollaron en primer lugar antes de salir al mercado civil. De hecho, muchas veces, a la hora de indagar más a fondo sobre el funcionamiento de algunos sensores aplicados en estos vehículos, no se disponía de la información suficiente para conocer en detalle las características del funcionamiento del sistema, debido a la confidencialidad tanto de empresas como de militares. En este caso de los sensores, son fundamentales los sistemas de navegación como pueden ser la IMU y el GPS. Se puede decir que el resto de sensores son más bien opcionales, los cuales se usarán cuando se diseñe un dron específico para unos fines específicos. Pero los sistemas de navegación tanto inerciales como por satélites son, como ya se ha dicho anteriormente, un “must” dentro del mundo de los UAVs. Es verdad que la tecnología avanza a pasos agigantados, pero en este estudio se ha comprobado que existen dispositivos muy utilizados a día de hoy con problemas de precisión. Por tanto, es necesario seguir trabajando para conseguir reducir aún más sus fuentes de error para obtener datos que sean lo más precisos y seguros posible. No solamente se necesita desarrollar en particular la estructura y/o diseño de un sistema, sino que se deben desarrollar métodos de integración de sistemas que sean más fiables, capaces de minimizar los errores y que se aproximen lo más posible a la realidad que se desea medir. La actual implementación de nuevos sensores y sistemas en distintos dispositivos gracias a los avances tecnológicos hace necesario que el estudio y desarrollo de la integración de señales no se quede atrás. Todavía existen numerosos aspectos en los que el filtro de Kalman puede ser mejorado, es por ello que a día de hoy se siguen publicando estudios en los cuales se proponen nuevos métodos de implementación de dicho filtro. Es necesario conocer correctamente los principios en los que se fundamenta el algoritmo utilizado por este filtro, para así vislumbrar los posibles errores que se comete al realizar el proceso iterativo y así, poder mejorarlo desarrollando una técnica alternativa. Se ha observado también la facilidad con la que se encuentra actualmente una amplia gama de componentes electrónicos, desde baratos que te ofrecen un servicio básico hasta los más caros que cuentan con unos rangos de errores bastante reducidos y abarcan un mayor espectro dentro de las actuaciones que puede llevar a cabo. De hecho, mismo en los smartphones que se utilizan cada día, se encuentran por ejemplo el GPS, acelerómetros, giróscopos y magnetómetros. Con esto se pretende señalar que los sistemas de navegación no son algo que se encuentre fuera de nuestro alcance, sino que 76 Conclusiones y Líneas Futuras además estamos en contacto con ellos a diario. En cuanto al caso práctico llevado a cabo en este estudio del Trabajo de Fin de Grado, se han encontrado dificultades al intentar materializar los conceptos teóricos. La literatura que se ha utilizado para entender los algoritmos en los que se basa el filtro de Kalman es bastante confusa. Ha sido todo un reto entender la información encontrada, descubriendo que para el desarrollo de problemas parecidos, existen diversas formas de interpretar el filtro de Kalman. Se recomienda que en estudios posteriores, se detalle en profundidad la funcionalidad de cada parámetro utilizado, y por qué adopta la forma que tiene. También se ha observado claramente que el GPS no proporciona unos datos demasiado fiables en cuanto a la variable altitud, tal y como se puede apreciar en las gráficas representadas en el apartado del Caso Práctico. La dependencia que presenta este sistema de una organización estadounidense, anima al resto de la comunidad tecnológica a desarrollar sus propios métodos de navegación basados en GNSS, como son la europea con Galileo, la rusa con GLONASS, la china con BeiDou, etc. La integración de varios sistemas de este tipo aumentaría notablemente la calidad de la señal recibida, y por tanto, la fiabilidad de la información, puesto que a día de hoy se encuentra muy sujeta a las decisiones que tomen las autoridades estadounidenses sobre este ámbito. Como líneas futuras para este trabajo, sería interesante implementar el filtro de Kalman para señales procedentes no solamente del GPS, sino también de una IMU. Como se dijo en el apartado del estudio sobre el GPS, es recomendable, cuando un vehículo se encuentra navegando, disponer de sistemas de medidas inerciales aparte de sistemas GNSS. De hecho, la IMU es el sistema inercial más extendido dentro del ámbito de los UAVs. Típicamente los drones que se encuentran en el mercado cuentan con GPS+IMU, por lo tanto, aplicar un filtro de Kalman a las señales recibidas por dichos sistemas sería un problema que despierta bastante interés por el extendido uso que se le podría dar. De esta forma, sería posible estudiar el caso en el que ambas señales son obtenidas de sensores reales e integradas. Es decir, la señal con la que se ha trabajado en este estudio obtenida mediante Google Maps sería sustituida por una real. En ese caso, la simplificación de trayectoria rectilínea y uniforme llevada a cabo en este documento no sería necesaria. En el presente trabajo, se ha simulado en Matlab un filtro de Kalman sencillo para esclarecer los conceptos que se utilizan en dicho algoritmo, pero existen otros tipos de filtros también utilizados en la literatura relacionada. Por ejemplo, [6] propone utilizar un método denominado “Square Root Cubature Kalman Filter, CKF” con el fin de reducir los problemas relacionados con el filtro de Kalman usual: cantidad de cálculos a realizar, complejidad de las fórmulas y las transformaciones, “baja” exactitud, poca convergencia o incluso divergencia. Mediante su estudio, se demuestra que utilizando este método se obtienen mayores estabilidad y precisión que utilizando el filtro de Kalman y otra modalidad llamada “Unscented Kalman Filter, UKF”. También según [6], este último UKF atrae la atención en este área puesto que obtiene una media y una covarianza usando un muestreado determinista, pero como inconveniente pero tarda más tiempo que por ejemplo el “Extended Kalman Filter, EKF” [23]. Además, cuando el número de estados es mayor de 3, el UKF no funciona correctamente, incluso puede llegar a detener su operación. Se propone estudiar alguna de estas variantes del filtro de Kalman para implementarla en un sistema que integre señales IMU y GPS. Como se ha recalcado varias veces durante este trabajo, la navegación de cualquier tipo de vehículo es normalmente un problema de estimación de estado altamente no lineal. El EKF ha jugado un papel importante en dicha estimación durante décadas [1]. Sería interesante intentar implementar estas versiones del filtro de Kalman y comparar los resultados con los obtenidos con otros métodos de filtrado. 77 Apéndice A Varianza & Covarianza Medida individual: xi Media de las medidas: ¯x Desviación de la media: ¯x−xi Cuadrado de la desviación: (¯x−xi)2 Varianza: σ2 x=PN i=1(¯x−xi)2 N Covarianza: σxσy=PN i=1(¯x−xi)(¯y−yi) Desviación Estándar: σx=pσ2 x=qPN i=1(¯x−xi)2 N 2D=σ2 xσxσy σyσxσ2 y 3D=  σ2 xσxσyσxσz σyσxσ2 yσyσz σzσxσzσyσ2 z   78 Apéndice B Código Matlab de caso práctico Implementación del filtro de Kalman % Jose Luis Carretero Rodriguez % Trabajo Fin de Grado - version 3D % S e al GPS/Google Maps con Filtro de Kalman % Nomenclatura: % LLH = Latitud, Longitud, Altitud % [x;y;vect_alt] = datos GPS % [x_est;y_est;altitud_est] = datos GoogleMaps % [vect_x_kalman;vect_y_kalman;vect_alt_kalman] = datos filtro clear; clc; format long %% Iniciar variables cont=0; any_south_lat=0; any_west_long=0; vector_lat=[]; vector_long=[]; vect_alt=[]; x=[];y=[];z=[]; velocidad=[]; vector_kalman=[]; % Modelo de la Tierra WGS-84 a=6378137; f=1/298.257223563; %% Obtengo datos del fichero NMEA TLINE = fopen('puente1.log'); tline = fgets(TLINE); while ischar(tline); token=strtok(tline); % VELOCIDADES % Obtengo la velocidad a partir de la linea GPVTG if strcmp(tline(2:6),'GPVTG')==1 [~,~,~,~,~,~,speed]=NMEA_VTG(tline); % Compruebo que no hay errores en el mensaje flag=checksum(token); if flag==0 disp('Error')% Si hay error, muestro por pantalla else % Convierto la cadena a n mero speed=str2num(speed); % Almaceno las velocidades velocidad=[velocidad; speed]; % Si no hay error, guardo end % LATITUD LONGITUD ALTITUD % Obtengo la posicion a partir de la linea GPGGA elseif strcmp(tline(2:6),'GPGGA')==1 flag=checksum(token); if flag==0 disp('Error') end [utc,lat,NS,long,EW,~,~,~,alt,unit_alt]=NMEA_GGA(tline); alt; if strcmp(NS,'S')% Si estoy en el Sur, mostrar %disp('Hay sur') % Nunca estoy en el Sur en este caso lat=num2str(-str2double(lat)); any_south_lat=1; end if strcmp(EW,'W')% Si estoy en el Oeste, mostrar %disp('Hay Oeste') % Siempre estoy en el Oeste en este caso long=num2str(-str2double(long)); % Cambiar signo any_west_long=1; 79 B Código Matlab de caso práctico end % Guardo todas las medidas de posicion. LLH % Formato string vector_lat=[vector_lat; lat]; if size(long)~=9 long=strcat(long,'0'); % dimensiones correctas end vector_long=[vector_long; long]; % Formato numero vect_alt=[vect_alt; str2double(alt)]; end tline=fgets(TLINE); % Obtener nueva l nea end % Convertir lat y long de string a numero for i=1:length(vector_lat) vect_lat(i,1)=str2double(vector_lat(i,1:2))+... str2double(vector_lat(i,3:end))/60; vect_long(i,1)=str2double(vector_long(i,2))+... str2double(vector_long(i,3:end))/60; end % Convertir LLH->xyz for i=1:length(vector_lat) xyz=llh2xyz([vect_lat(i) -vect_long(i) vect_alt(i)],a,f); % Datos de posicin del GPS (x,y,z) x=[x; xyz(1)]; y=[y; xyz(2)]; z=[z; xyz(3)]; end % Calcular altitud media alt_media=sum(vect_alt)/66; % plot(x) % plot(y) % plot(vect_alt); hold on % plot(alt_media*ones(66,1)) % plot3(x,y,z,'LineWidth',2),grid % xlabel('Eje X') % ylabel('Eje Y') % zlabel('Eje Z') % Matriz de covarianza R - error en GPS % Primero calculo la media del error z sumatorio_GPS=[0;0;0]; for i=1:65 sumatorio_GPS=sumatorio_GPS+... [x(i+1)-x(i);... y(i+1)-y(i);... vect_alt(i+1)-vect_alt(i);]; end media_z=sumatorio_GPS/66; % Calculo R R1=0*ones(3); for i=1:65 R1=R1+[x(i+1)-x(i)-media_z(1);... y(i+1)-y(i)-media_z(2);... vect_alt(i+1)-vect_alt(i)-media_z(3);]*... [x(i+1)-x(i)-media_z(1);... y(i+1)-y(i)-media_z(2);... vect_alt(i+1)-vect_alt(i)-media_z(3);]'; end R=R1/65; % 65 segundos dura la trayectoria - 65 muestras. 80 B Código Matlab de caso práctico %% Google Maps % Puente segun Google Maps (lat,long,alt) lat_ini_est=38.878408; long_ini_est=-6.980620; lat_fin_est=38.883509; long_fin_est=-6.984960; alt_est=167; % Convierto a (x,y,z) est_ini=llh2xyz([lat_ini_est long_ini_est alt_est],a,f); x_ini_est=est_ini(1); y_ini_est=est_ini(2); est_fin=llh2xyz([lat_fin_est long_fin_est alt_est],a,f); x_fin_est=est_fin(1); y_fin_est=est_fin(2); % Calculo velocidades - (final-inicial)/tiempo velocidad_x=(x_fin_est-x_ini_est)/65; velocidad_y=(y_fin_est-y_ini_est)/65; % Creo vector de componentes de posicion for t=1:66 x_est(t)=x_ini_est+velocidad_x*t; y_est(t)=y_ini_est+velocidad_y*t; altitud_est(t)=alt_est+0*t; end x_est=x_est'; y_est=y_est'; altitud_est=altitud_est'; % Matriz de covarianza Q % Primero el ruido medio de w sumatorio_est=[0;0;0]; for i=1:65 sumatorio_est=sumatorio_est+... [x_est(i+1)-x_est(i);... y_est(i+1)-y_est(i);... altitud_est(i+1)-altitud_est(i);]; end media_w=sumatorio_est/66; % Matriz Q Q1=0*ones(3); for i=1:65 Q1=Q1+... [x_est(i+1)-x_est(i)-media_w(1);... y_est(i+1)-y_est(i)-media_w(2);... altitud_est(i+1)-altitud_est(i)-media_w(3);]*... [x_est(i+1)-x_est(i)-media_w(1);... y_est(i+1)-y_est(i)-media_w(2);... altitud_est(i+1)-altitud_est(i)-media_w(3);]'; end Q=Q1/65; % figure('Name','Posicin x') % plot(x_est,'r'); hold on % plot(x,'b') % % figure % plot(y_est,'r'); hold on % plot(y,'b') % % figure % plot(x,y,'r'); hold on % plot(x_est,y_est,'b') % xlabel('Eje X') % ylabel('Eje Y') 81