Control robusto no lineal de un satélite en presencia de perturbaciones e incertidumbre en los parámetros
Abstract
En este proyecto se va a simular leyes de control de la orientación de un satélite en una órbita baja. La dinámica del satélite se va a modelar con las ecuaciones clásicas de rotación y va a estar influída por el efecto de la gravedad de la Tierra, la cual crea una perturbación sobre el sistema. El desarrollo por tanto de este proyecto va a tener varias partes. Se explicará la tecnología que hace que esto sea posible, la teoría de control sobre la que se basa la simulación, se explicará también las leyes de control. Este proyecto se basa en lo desarrollado en el libro “Robust Autonomous Guidance”.
Full text
Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Control robusto no lineal de un satélite en presencia de perturbaciones e incertidumbre en los parámetros Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Miguel Díaz Cabrera Tutor: Dr. Eduardo Fernández Camacho Sevilla, 2019
iii Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Control robusto no lineal de un satélite en presencia de perturbaciones e incertidumbre en los parámetros Autor: Miguel Díaz Cabrera Tutor: Dr. Eduardo Fernández Camacho Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019
v Trabajo de fin de grado: Control robusto no lineal de un satélite en presencia de perturbaciones e incertidumbre en los parámetros Autor: Miguel Díaz Cabrera Tutor: Dr. Eduardo Fernández Camacho El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2019 El Secretario del Tribunal
Contents 1Introducción 4 1.1 Historia .............................. 4 1.1.1 Sputnik1 ......................... 4 1.1.2 Explorer1......................... 5 1.2 Sistemas de un satélite . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Clasificación de las Órbitas . . . . . . . . . . . . . . . . . . . . 6 1.4 Clasificación de satélites por misión . . . . . . . . . . . . . . . 7 1.5 Industria.............................. 8 1.5.1 Producción ........................ 9 2Tecnologíadesensoresyactuadoresdeunsatélite 10 2.1 Sensores.............................. 11 2.1.1 Giroscopios ........................ 11 2.1.2 Sensores solares . . . . . . . . . . . . . . . . . . . . . . 12 2.1.3 Detectores de estrella “Star tracker” . . . . . . . . . . . 12 2.1.4 Magnetómetro . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Actuadores ............................ 12 2.2.1 Propulsores ........................ 12 2.2.2 Ruedas de reacción . . . . . . . . . . . . . . . . . . . . 12 2.2.3 Velassolares........................ 13 3Conceptos 13 3.1 Posición y velocidad . . . . . . . . . . . . . . . . . . . . . . . 13 3.1.1 Matriz de rotación . . . . . . . . . . . . . . . . . . . . 13 3.1.2 Cuaterniones . . . . . . . . . . . . . . . . . . . . . . . 14 3.1.3 Cinemática de la Orientación: Matriz de rotación y cuaterniones........................ 15 3.2 Conceptos matemáticos . . . . . . . . . . . . . . . . . . . . . . 16 3.3 Sistema Lineal y No Lineal . . . . . . . . . . . . . . . . . . . . 16 3.4 Inmersión de un sistema . . . . . . . . . . . . . . . . . . . . . 16 3.5 Estabilidad ............................ 17 3.5.1 Puntos de equilibrio . . . . . . . . . . . . . . . . . . . 17 3.5.2 Estabilidad de puntos de equilibrio . . . . . . . . . . . 17 3.5.3 Teorema de Lyapunov . . . . . . . . . . . . . . . . . . 18 3.5.4 Teorema de LaSalle . . . . . . . . . . . . . . . . . . . . 19 4 Teoría del problema 20 4.1 Caso en que se tiene información plena: conocimiento del estado x(t)ydelaperturbaciónw(t)............... 21 1
4.2 Caso en que solo se tiene información parcial: conocimiento de la variable de error e(t).................... 21 4.2.1 Problema de control . . . . . . . . . . . . . . . . . . . 21 4.2.2 Solución al problema de control . . . . . . . . . . . . . 22 4.2.3 Ecuaciones de regulación . . . . . . . . . . . . . . . . . 22 4.2.4 Propiedad de Modelo Interno . . . . . . . . . . . . . . 23 4.2.5 Interpretación geométrica de las ecuaciones de regulación 23 5Sistemadereferencia 24 5.1 Centrado en la Tierra . . . . . . . . . . . . . . . . . . . . . . 24 5.2 Centrado en el satélite . . . . . . . . . . . . . . . . . . . . . . 25 5.3 ÓrbitaLocal ........................... 25 5.4 Posición y velocidad angular respecto al sistema de coordenadas 26 6Perturbacionesqueafectanalsistema 27 6.1 Gradiente gravitacional . . . . . . . . . . . . . . . . . . . . . 27 6.1.1 Fuerza gravitacional . . . . . . . . . . . . . . . . . . . 27 6.1.2 Torque gradiente gravitacional . . . . . . . . . . . . . . 28 7Dinámicadelsistema 29 7.1 Física del satélite . . . . . . . . . . . . . . . . . . . . . . . . . 29 7.2 Ecuaciones dinámicas del satélite . . . . . . . . . . . . . . . . 31 7.3 Ecuaciones dinámicas de error del satélite . . . . . . . . . . . 31 8LeydecontrolNolineal 34 8.1 Acción por Modelo Interno uim ................. 34 8.2 Acción Estabilizadora ust .................... 35 8.2.1 Estabilidad de la interconexión del error de la dinámica yelmodelointerno . . . . . . . . . . . . . . . . . . . . 36 8.3 Ley de control: suma de la acción por modelo interno uim y acción estabilizadora ust ..................... 36 9SimulaciónusandoSimulink 37 9.1 Misión: Vuelo inercial . . . . . . . . . . . . . . . . . . . . . . . 37 9.1.1 Objetivo de la misión . . . . . . . . . . . . . . . . . . . 37 9.1.2 Diseño de la acción de control por modelo interno uim .41 9.1.3 Estabilidad de la interconexión . . . . . . . . . . . . . 44 9.1.4 Simulación del algoritmo de control en Simulink . . . . 45 9.2 Misión: Vuelo orientado hacia la tierra . . . . . . . . . . . . . 51 9.2.1 Objetivo de la misión . . . . . . . . . . . . . . . . . . . 51 9.2.2 Diseño de la acción de control por modelo interno uim .52 2
9.2.3 Simulación del algoritmo de control . . . . . . . . . . . 53 9.3 Conclusión............................. 56 3
1Introducción En este proyecto se va a simular leyes de control de la orientación de un satélite en una órbita baja. La dinámica del satélite se va a modelar con las ecuaciones clásicas de rotación y va a estar influída por el efecto de la gravedad de la Tierra, la cual crea una perturbación sobre el sistema. El desarrollo por tanto de este proyecto va a tener varias partes. Se explicará la tecnología que hace que esto sea posible, la teoría de control sobre la que se basa la simulación, se explicará también las leyes de control. Este proyecto se basa en lo desarrollado en el libro “Robust Autonomous Guidance”. Figure 1: Satélite 1.1 Historia El primer satélite en entrar en servicio fue el artefacto aeroespacial Sputnik 1delaUniónSoviéticael4octubrede1957iniciandoelprogramasoviético Sputnik. Esto fue el comienzo de la carrera espacial entre Estados Unidos y la Unión Soviética. 1.1.1 Sputnik 1 El Sputnik 1 fue un satélite que fue lanzado en una órbita elíptica baja y que estuvo en servicio durante tres semanas hasta que las baterías se agotaron. Este satélite tenía forma de bola metáica con un diámetro de 58 cm, con 4antenasderadioparaemitirseñales. Sirvióparaestimarladensidadde la capa alta de la atmósfera gracias a la fuerza de arrastre del aire en su órbita. Adicionalmente la propagación de las ondas de radio proporcionó información sobre la ionosfera. 4
Figure 5: Space shuttle siempre orientada hacia la Tierra, o que el calor y frío aportado por la radiación solar se aproveche de una manera concreta, o por guiado: pequeñas maniobras a propulsión se deben ejecutar en la correcta dirección. 2.1 Sensores Los sensores generan información que refleja el ritmo de cambio en la orientación. Estos tienen siempre un umbral de medida. 2.1.1 Giroscopios Son dispositivos que miden la rotación en un espacio tridimensional sin tener que tener una referencia. Para ello se basan en el principio de conservación del momento angular. Tradicionalmetne, un giroscopio consiste en una rueda montada sobre dos o tres suspensiones cardán Hoy en día existen sistemas electromecánicos (MEMS) que son giroscopios minituarizados en dispositivos electrónicos. 11
2.1.2 Sensores solares Es un instrumento de navegación que detecta la posición del Sol. Miden el ángulo respecto al Sol y lo indican mediante señales continuas y discretas, respectivamente. 2.1.3 Detectores de estrella “Star tracker” Es un dispositivo óptico que mide la posición de las estrellas utilizando cámras. Teniendo en cuenta que ya se conocen las posiciones de muchas estrellas con gran precisión este dispositivo montado sobre un satélite se puede utilizar para conocer la orientación respecto a las estrellas. Para esto el detector de estrella tiene que obtener una imagen de las estrellas, medir la posición respecto al sistema referencial del satélite e identificar las estrellas de manera que un procesador pueda comparar con la bases de datos la posición absoluta de estas estrellas. 2.1.4 Magnetómetro Es un instrumento que mide el campo magnético la dirección de este. Luego se compara con un mapa del campo magnético terrestre guardado en una memoria a bordo. Y si sabe la posición del satélite se puede medir la orientación. 2.2 Actuadores 2.2.1 Propulsores Un propulsor es capaz de crear un momento respecto al centro de gravedad del satélite para cambiar la orientación girando sobre los 3 ejes principales. Para esto se deben distribuir de manera organizada en la estructura del satélite de manera que puedan generar momentos de fuerza sin generar una translación del vehículo. Las limitaciones de estos son el uso de combustible, gasto del motor, y de los ciclos de las válvualas. La eficiencia del combustible se mide con el menor impulso que pueda aportar el propulsor. 2.2.2 Ruedas de reacción Estos son ruedas que giran debido a un motor eléctrico equipado en el satélite. Cuando estas ruedas se activan hacen por el conservación del momento angular el satélite empiece a girar en el sentido opuesto. 12
Debido a su baja masa y a que son controladas por ordenador tienen una gran precisión lo que permite realizar pequeños giros. 2.2.3 Velas solares Son dispositivos que producen una fuerza debido a luz que impacta contra estas. Este tipo de dispositivos ahorra combustible en una misión de larga duración debido a que no tiene un gasto energético. 3 Conceptos 3.1 Posición y velocidad 3.1.1 Matriz de rotación Un sistema de coordenadas Fse define como F={O, !i, !j, !k}donde se define un centro de coordenadas y una base canónica en R3, i=0 @ 1 0 01 A,j=0 @ 0 1 01 A,k=0 @ 0 0 11 A Una matriz de rotación expresa la diferencia de posición angular u orientación entre dos sistemas de coordenadas distintas. Se llama Fa={Oa,!ia,!ja,! ka} aunsistemayFb={Ob,!ib,!jb,!kb}al otro sistema. La orientación de Fb relativa a a Faviene expresada por la matriz de rotación: Rab =0 B @ !ib·!ia !jb·!ia !kb·!ia !ib·!ja !jb·!ja !kb·!ja !ib·! ka !jb·! ka !kb·! ka 1 C A Esta matriz de rotación es un elemento del grupo especial ortogonal SO(3) ⇢R3x3,esdecir,elconjunto: SO(3) = {R2R3x3:RRT=I,det(R)=1} Yrelacionadosvectoresunitariosexpresadosencadaunodeestosdos sistemas de coordenadas: ⇣!ib !jb !kb⌘=Rab ⇣!ia !ja ! ka⌘ Un vector genérico !vse resuelve en FayFbcomo 13
!v=va 1 !ia+va 2 !ja+va 3 ! ka y !v=vb 1 !ib+vb 2 !jb+vb 3 !kb La relación entre los dos vectores: va=0 @ va 1 va 2 va 31 A,v b=0 @ vb 1 vb 2 vb 31 A viene dada por la matriz de rotación: va=Rabvb 3.1.2 Cuaterniones Otra manera de caracterizar la orientación entre dos sistemas de coordenadas viene dado por el uso de cuaterniones. Estos son un cuadruple de números reales (q0,q 1,q 2,q 3)que cumple con la restricción: Xq2 i=1 Yquevienedefinidoporelconjunto: S4={x2R4:nxn=1} Normalmente cada cuaternión se expresa con un número escalar y una parte vectorial de dimensión 3: q=✓q0 q◆ donde q0es la parte escalar y q=0 @ q1 q2 q31 A es la parte vectorial. La matriz de rotación definida anteriormente es sencilla de visualizar, sin embargo, los cuaterniones no son facilmente visualizables ya que se habla de una dimensión 4. Lo que se debe saber es que cada conjunto de 4 números 14
representa una orientación particular, son fáciles de manipular y de tratar computacionalmente. Para cada cuaternión qexiste una matriz concreta expresada por: R(q)=0 @ 12q2 22q2 32q1q22q0q32q1q3+2q 0q2 2q1q2+2q 0q312q2 12q2 32q2q32q0q1 2q1q32q0q22q2q3+2q 0q112q2 12q2 21 A que satisface RT(q)R(q)=Iydet(R(q)) = 1.De igual manera s puede demostrar que para cada matriz de rotación Rexiste un cuaternión qde manera que: R=R(q) 3.1.3 Cinemática de la Orientación: Matriz de rotación y cuaterniones Si un sistema de coordenadas Fbvaría en el tiempo respecto a Fala matriz de rotación Rab varía en el tiempo. Esta relación viene dada por ˙ Rab =RabSkew(wb ab)=Rab ⇥wb ab Igualmente usando cuaterniones esta relación se puede expresar como ˙q =1 2E(q)w donde E(q)=✓qT q0I+Skew(q)◆=0 B B @ q1q2q3 q0q3q2 q3q0q1 q2q1q0 1 C C A ydondeseintroduceeloperadorlinealSkew() por conveniencia. De manera que Skew(v)=0 @ 0v3v2 v30v1 v2v101 A Este operador es equivalente a decir ”⇥v” 15
3.2 Conceptos matemáticos •Una proposición es una declaración que es verdadera o falsa. Es similar aunteoremaaunqueelteoremasesueleutilizarmáscomounresultado principal. •Una declaración que es de interés para demostrar un teorema se llama lemma. 3.3 Sistema Lineal y No Lineal Un sistema lineal viene representado por ˙x=Ax +Bu y=Cx +Du Un sistema no lineal viene representado por: ˙x=f(x, u) y=k(x, u) 3.4 Inmersión de un sistema Definición: Dado dos sistemas con la misma salida y: (˙x=f(x),x2 y=h(x),y 2Rm y (˙ X=F(X),X 2X Y=H(X),Y 2Rm se dice que {,f,h}está inmerso dentro de {X,F,H}si existe una función ⌧:!Xque cumple con ⌧(0) = 0 y @⌧ @xf(x)=F(⌧(x)) h(x)=H(⌧(x)) para todo x2 16
La razón por la que este concepto es importante para el problema de regulación que se resuelve en este proyecto es porque da la posibilidad de tener el sistema autónomo: ˙w=s(w) u=c(w) inmerso dentro de otro sistema: ˙ ⇠='(⇠) u=(⇠) que tenga una propiedades deseadas. 3.5 Estabilidad 3.5.1 Puntos de equilibrio Un punto x⇤llama punto de equilibrio si tiene la propiedad que para cualquier estado del sistema x,concondicióninicialx(0) = x⇤permanece en x⇤para todo el t. 3.5.2 Estabilidad de puntos de equilibrio Aparte del concepto del punto de equilibrio interesa saber el concepto de estabilidad de ese punto de equilibrio. De manera general la estabilidad del punto de equilibrio implica que si el estado xse encuentra cerca del punto de equilibrio x⇤,permanecerácercaaesteparatodoeltiempofuturot>0. Formalmente se define: “Elpuntodeequilibriox=0es 1. estable si para cualquier ✏>0,se puede encontrar un >0de manera que kx(t0)k<implica que kx(t)k<✏para todo t t0 2. inestable si no es estable 3. asintóticamente estable si es estable y se puede encontrar un >0de manera que kx(t0)k<implica que lim t!1 x(t)=0 17
Figure 6: equilibria 4. globalmente asintóticamente estable si es estable y lim t!1 x(t)=0, para todo x(t0) ” La definición 3 y 4 refleja como se aplican estos conceptos a este proyecto. Si se define la orientación y la velocidad angular del satélite como el estado xyexistenperturbacionesquehacenqueesteestadosealejedelpuntode equilibrio definiendo esta distancia por se quiere que se vuelva al equilibrio. La definición 3 aplica al caso local, es decir para un punto cercano al equilibrio. La definición 4 es mucho más estricta y aplica para cualquier punto del entorno global. 3.5.3 Teorema de Lyapunov Se define el teorema de estabilidad de Lyapunov como: “Considerando el sistema ˙x=f(x) donde x2Rnyasumiendoqueexistesoluciónparax(t)en un conjunto abierto D⇢Rnque contiene el origen (esto es 02D).Sea x=0un punto 18
de equilibrio del sistema. Y considerando una función continua diferenciable V(x):D!Rde manera que V(0) = 0 yV(x)>0para x2Dcon x6=0. Suponiendo ahora que para la solución x(t), ˙ V(x)0 Entonces, x=0es estable. De hecho, si ˙ V(x)0para x 6=0, en las soluciones x(t),entoncesx=0es asintóticamente estable.” Adicionalmente se puede presentar otro resultado similar para demostrar estabilidad global asintótica: “SuponiendoqueelteoremaanteriorsecumpleyquelafuncióndeLyapunov V(x)tiene la propiedad de que kxk!1 implica que V (x)!1 Luego x=0es globalmente asintóticamente estable.” 3.5.4 Teorema de LaSalle “Suponiendo un sistema ˙x=f(x) donde x2Rn.Se dice que un conjunto M⇢Rnes invariante respecto a este sistema si x(0) 2M hace que x(t)2M para todo t 2R ” Esto significa que si la solución x(t)permanece a Mpara un instante de tiempo, sigue permaneciendo en Mpara todo t. Interesa saber la idea de lo que se quiere expresar cuando se dice que un punto se acerca a un conjunto: “Sedicequex(t)se acerca un conjunto Mcuando t!1si para cada ✏>0existe un T>0tal que dist(x(t),M)<✏, para todo t T, 19
donde dist(p, M)denota la distancia de un punto paunconjuntoM, que se define como dist(p, M)=min kpqk ” El Teorema de LaSalle es similar al Teorema de Lyapunov y permite establecer argumentos de estabilidad en conjuntos: “Considerando el sistema ˙x=f(x) yasumiendoqueexisteunasoluciónx(t)alaecuacióndeestadosenun conjunto abierto D⇢Rn.Suponiendo que para una solución x(t)se puede encontrar un conjunto cerrado y acotado ⌦⇢Dde manera que x(t)2⌦ para todo t0.SeaunaV(x):D!Runa función continua diferenciable de manera que a lo largo de las trayectorias x(t)la derivada ˙ V0en ⌦.Sea Eel conjunto de puntos dentro de ⌦donde ˙ V(x)=0,esto es E={x2⌦: ˙ V(x)=0}.Sea Mel conjunto invariante más grande del sistema contenido en E. Entonces x(t)!Mcuando t!1” 4 Teoría del problema Dado el sistema: ˙x=f(x, w, u) e=h(x, w) El problema de control que se resuelve en este proyecto es un problema clásico de control de regulación de la trayectoria. Es decir, conseguir a través de una retroalimentación que la salida del sistema y(t)coincida con una referencia deseada yref (t)al igual que rechazar, si existe, una perturbación no deseada w(t). En cualquier caso, el problema se expresa en conseguir que el error de regulación e(t)=yref (t)y(t) llegue a 0. 20
Figure 10: Atracción de la gravedad 6Perturbacionesqueafectanalsistema Acontinuaciónsevequeformatieneloquesehallamadowyqueintroduce en el sistema unas perturbaciones no deseadas. 6.1 Gradiente gravitacional 6.1.1 Fuerza gravitacional El torque por el gradiente gravitacional se debe al hecho que la fuerza de gravedad debido a la masa de la Tierra ejercida sobre el satélite no es constante con la distancia de la distancia sino que crece cuadraticamente. En consecuencia, la parte del satélite que se encuentra más lejos de la tierra se ve menos afectada por la Tierra que la parte más cercana a esta. Esta diferencia de fuerzas crea un momento respecto al centro de masas. La fuerza gravitacional sobre un elemento diferencial de masa dm localizado a una distancia !⇢del centro de masas del satélite viene dada por d!f=µ(!r+!⇢) |!r+!⇢|3dm donde !res la posición orbital del centro de masas del satélite. Vamos a suponer que la distancia ⇢es mucho menor que la distancia res decir ⇢⌧r 27
esto es que la mayor distancia de un punto del satélite al centro de masas es mucho menor que la distancia de la Tierra al centro de masas del satélite. Y donde µ=GM. Por tanto podemos aproximar utilizando la expansión de las series de Taylor el término: |!r+!⇢|3⇡1 r3(1 3!r·!⇢ r2) Sustituyendo está expresión en la fuerza que antes hemos indicado: d!f=µ(!r+!⇢) r3(1 3!r·!⇢ r2)dm 6.1.2 Torque gradiente gravitacional El torque total sobre el centro de masas viene dado por: !⌧grav =3µ r5ZV !⇢⇥d!f= =ZV !⇢⇥(µ(!r+!⇢) r3)dm +ZV !⇢⇥(3µ!r·!⇢ r5(!r+!⇢))dm donde la integral se calcula sobre el volumen entero del satélite. Teniendo en cuenta RV !⇢dm =0ya que el origen está en el centro de masas. Y que !⇢⇥!⇢=0,laexpresióndeltorquesesimplificaa: !⌧grav =3µ r5ZV !⇢⇥!r(!r·!⇢)dm Debido a la distribución de masas del satélite la fuerza de la gravedad que actúa sobre el satélite ejerce mayor fuerza en la parte del satélite más cercana al satélite, esto crea un torque que pueda afectar a la posición deseada del satélite. Se puede llegar a la conclusión que este momento viene dado por: !⌧grav =3µ r3(!k0⇥!J·!k0) donde !k0es el vector unitario de la base de órbita local. Si expresamos este torque respecto al sistema de coordenadas fijadas al centro de masas del satélite Fb,este vector unitario se puede expresar: !k0=Rbo !kb=RT ob !kb 28
donde Rob es la matriz de orientación que relacióna la base F0con Fb. Por tanto este torque respecto a este sistema de coordenadas mencionado quedaría como: !⌧b grav =3µ r3(RT ob !k)⇥!J·RT ob !k 7Dinámicadelsistema 7.1 Física del satélite El satélite que se considera en este proyecto está en una órbita circular afectado por la gravedad de la Tierra y se estudia su movimiento mediante las leyes de Newton que relacionan el movimiento de una partícula. La 1ºLey dice que si no hay ninguna fuerza actuando sobre un cuerpo, el cuerpo en reposo permanecerá en reposo y el cuerpo en movimiento permancerá en movimiento constante en una línea recta. La 2ºLey dice que si se aplica una fuerza habrá un cambio de velocidad, es decir, una aceleración proporcional a la magnitud de la fuerza y en la dirección que se aplica la fuerza. Esta ley se expresa como : F=ma donde Fes la fuerza, mes la masa de la partícula, y aes la aceleración. La 3ºLey dice que si un cuerpo provoca una fuerza sobre otro cuerpo, este otro cuerpo crea otra fuerza de igual magnitud pero dirección opuesta sobre el cuerpo inicial. En la Ley de gravitación universal Netwon expresa que dos partículas con masas m1ym2yseparadasporunadistanciarestán atraídas por cada una con igual fuerza y en dirección a la línea que va de una a la otra. Esta magnitud de la fuerza es: F=Gm1m2 r2 donde G=6.67259 ⇥1011Nm2 kg2es la constante de gravitación. En una órbita circular como es el caso el satélite se encuentra en caída libre. El satélite acelera hacia el centro de la Tierra mientras se mueve en una línea recta. La velocidad de este cambia continuamente en dirección pero no en magnitud. A partir de las Leyes de Netwon se ve que la dirección de la velocidad está cambiando, existe aceleración. Esta aceleración, llamada aceleración centrípeta tiene dirección al centro de la Tierra y viene dada por 29
a=v2 r donde ves la velocidad del satélite y res el radio del círculo imaginario sobre el que gira. Esta aceleración seegún las Leyes de Newton debe ser creada por una fuerza, la fuerza centrípeta cuya magnitud viene dada por: F=mv2 r La dirección de Fen cualquier instante tiene igual dirección que a,esto es, hacia radialmente hacia dentro. Ahora haciendo balances de fuerzas sobre el satélite según la 2ºLey de Newton: fuerza centrípeta y fuerza de gravedad, se puede calcular la velocidad del satélite en equilibrio: v2 r=GM r2 o v=rGM r Ylavelocidadangulardelsatélitecomo: w=v r!w2=v2 r2 Teniendo en cuenta la expresión de la aceleración centrípeta anterior: w2=a r Yportanto, w2=GM r3 Se tiene el módulo de !web que es la velocidad con la que el satélite girar alrededor de la Tierra en función de la constante gravitacional, la masa de la Tierra y de la distancia de la Tierra respecto al satélite. 30
7.2 Ecuaciones dinámicas del satélite El movimiento que sigue el satélite vienen descritas por las ecuaciones de Euler para un sólido rígido en rotación, las cuales se pueden encontrar en cualquier texto estándar de mecánica clásica: ˙ Reb =Reb Skew(wb eb) J˙wb eb =wb eb ⇥Jwb eb +⌧b en donde la matriz de inercia es expresada en el sistema de coordenadas centradas en el satélite, es simétrica y positiva. Este modelo se ha derivado de las ecuaciones de movimiento de Euler para un sólido rígido. El término Reb es la posición relativa de los ejes de posición del cuerpo respecto al sistema referencial de la tierra. Viene expresada en términos matriciales, la podemos llamar matriz de rotación, y tiene varias propiedades que se serán útiles para el desarrollo del controlador. El segundo término wb eb es la velocidad angular del sátelite respecto a los ejes de referencia de la tierra. El operador Skew es el operador equivalente a Reb ⇥wb eb. Luego en la segunda ecuación, podemos observar el comportamiento dinámico del sólido implicando a las fuerzas que actúan sobre el sólido. El término J en este caso es la matriz de inercia del sátelite de dimensión 3 y que depende de la geometría y distribución de masas del mismo, está expresada respecto alosejesdelsatélite,loscualesvamosadefinircomolosejesprincipalesde inercia. Esta matriz Jserá uno de los párametros de diseño del sátelite. La ecuación dinámica depende también de la velocidad angular y del torque aplicado al cuerpo que viene expresado mediante ⌧. Este término !⌧=!⌧dist +!uestá compuesto de un término que agrupa los dos torques que van a causar perturbaciones en la dinámica del satélite, los cuales serán la gravedad, la cuál atrae más a la parte más “cercana” del satélite a su órbita, y el otro efecto serán el causado por el rozamiento con las partículas de aire presentes en la capa baja de la atmósfera que se han explicado con anterioridad. 7.3 Ecuaciones dinámicas de error del satélite Acontinuaciónsevaadeducirlasecuacionesdinámicasdelerrordelaposición del satélite. Para ello se va a tener en cuenta como se expresa el error de la orientación y el error de la velocidad angular. Se define a continuación primero, el error de la orientación como diferencia entre el la orientación real del satélite respecto a la órbita local o la Tierra y 31
de la orientación deseada del satélite respecto a la órbita local o la Tierra. Y segundo, se define el error de la velocidad angular como la diferencia entre la velocidad angular real que tiene el satélite respecto a la tierra y la velocidad angular deseada del satélite respecto a la Tierra. e R(t)=Rdb(t)=RT odRob =RT edReb ew=wb eb wb ed Como se ha dicho anteriormente se va a asumir que la ley de control tiene conocimiento sobre el error (e R, ew).Desarrollandolasecuacionesdinámicas se obtienen las ecuaciones dinámicas del error como: ˙ e R=e RSkew(ew) J˙ ew=⌧b 0+⌧b grav +ub donde se sustituye wb eb =ew+wb ed de manera que: ⌧b 0=ew⇥Jewwb ed ⇥Jewew⇥Jwb ed wb ed ⇥Jwb ed J˙wb ed La velocidad angular angular deseada del satélite respecto a la Tierra en el sistema de coordenadas del cuerpo del satélite es: wb ed =e RTwd ed Y en donde la velocidad angular deseada del satélite respecto a la Tierra en el sistema de coordenadas de la orientación deseada es la suma de: la velocidad angular deseada respecto a la órbita local y la de la velocidad angular de la órbita local respecto a la Tierra, ambas expresadas en el sistema de coordenadas de la orientación deseada: wd ed =wd od +wd eo =wd od +RT odwo eo yteniendoencuentaquelavelocidadangulardelaórbitalocalrespecto alatierraenelsistemadecoordenadasdelaórbitalocaleswo eo =w0j.Se tiene: wd ed =wd od +wd eo =wd od +RT odwo eo =wd od w0RT odj 32
Por tanto: wb ed =e RTwd ed =e RT[wd od w0RT odj] Yladerivadadewb ed se calcula utilizando las propiedades de la derivada del producto: ˙wb ed =d dt(e RT[wd od w0RT odj]) =ew⇥e RT[wd od w0RT odj]+ e RT[˙wd od w0˙ RT odj] Teniendo en cuenta que ˙wd od =0ya que se ha dicho que el satélite gira con velocidad constante: =ew⇥e RT[wd od w0RT odj]+w0e RT[wd od ⇥RT odj] En definitiva las ecuaciones dinámicas del error quedan expresadas: ˙ e R=e RSkew(ew) J˙ ew=⌧b 0(e R, ew, Rod,w od,µ)+⌧b grav(e R, Rod,µ)+ub en función de las variables de error y de los inputs exógenos. Por tanto si se consigue una acción de control ubtal que sea igual a: ub=⌧b 0(I,0,R od,w d od,µ)⌧b grav(I,Rod,µ) En el punto (e R, ew)=(I,0) de la dinámica del error existirá un punto de equilibrio del sistema: ˙ e R=e RSkew(ew) J˙ ew=0 En este punto se cumple que e R=Iyew=0para todo t>0.Paraellose debe conseguir crear un modelo interno que replique el comportamiento de este input deseado sin tener conocimiento del estado de las perturbaciones Rod,w od,µ. 33
8LeydecontrolNolineal 8.1 Acción por Modelo Interno uim El sistema para el cual se debe conseguir una regulación del error viene dado por: ˙ e R=e RSkew(ew) J˙ ew=⌧b 0+⌧b grav +ub donde el error en la orientación y la velocidad angular se hallan así. e R(t)=Rdb(t)=RT odRob ew=wob wod •La acción de control que cumple con las ecuaciones de regulación: @⇡ @ws(w)=f(⇡(w),c(w),w) 0=h(⇡(w),w) viene dada por ub=c(Rod,w d od,µ)=⌧b 0(I,0,R od,w d od,µ)⌧b grav(I,Rod,µ) c(Rod,w d od,µ)= (wd odw0RT odj)⇥J(wd odw0RT odj)+ w0J(wd od⇥RT odj)3w2 0(RT odk)⇥JRT odk haciendo que el punto (e R, ew)=(I,0) sea un equilibrio del sistema. En este equilibrio se define un mapa (e R, ew)=⇡(Rod,w od,µ)en el cual el estado (e R, ew)permanece a lo largo del tiempo cuando se aplica esta ley de control: ub=c(Rod,w d od,µ) yenelcuallaposiciónorbitalesequivalentealadeseadaRob =Rod al igual que la velocidad angular wob =wod 34
•El modelo interno que es capaz de recrear esta acción de control y cumple con las ecuaciones de modelo interno: @ @ws(w)=((w),w)) c(w)=✓((w),w)) es un sistema dinámico y se verá que forma tiene más adelante en la parte de simulación. Este sistema dinámico es capaz de reproducir la acción de control c(Rod,w d od,µ) necesaria que cumple con las ecuaciones de regulación sin tener conocimiento de los inputs exógenos (Rod,w d od,µ) el cual tiene una forma concreta dependiendo de la acción de control ub=c(Rod,w d od,µ)que pretenda reproducir, esto es dependiendo del tipo de misión. 8.2 Acción Estabilizadora ust La acción ust estabilizadora de control se diseña aplicando principios de estabilidad de Lyapunov y Teorema de La Salle. Primero se hace un cambio de variable: e ⇠=⇠(w) z=ew+k1eq donde k1es un parámetro de diseño. Luego se considera como función de Lyapunov según [1] a : V(e ⇠,eq, ew+k1eq)= 2e ⇠TPe ⇠+(1eq0)2+eqTeq+1 2(ew+k1eq)TJ(ew+k1eq) Ydespuésdevariasmanipulacionesmatemáticasyeligiendolaacción estabilizadora: ub st =k2(1 + n(ew+k1eq)n)( ew+k1eq) 35
Se llega la conclusión que la derivada de la función de Lyapunov esta superiormente acotada: ˙ V(e ⇠,eq,z) k1 2neq 2 n✏(1 + nzn)nz 2 n<0 donde ✏>0. Esta acción de control por argumentos de Lyapunov y LaSalle hace que el sistema en bucle cerrado sea estable. 8.2.1 Estabilidad de la interconexión del error de la dinámica y el modelo interno Después de crear el modelo interno ⇠que contrarresta el efecto del sistema exógeno y añadirlo como una variable más de estado se debe asegurar la estabilidad de la conexión (x, ⇠).Paraestoseañadeeltérminogst yquesediseña con argumentos de Lyapunov. Este término garantiza que la interconexión estado xymodelointerno⇠es estable. ˙ ⇠=⇠ +gst u=⇠ El valor de este término se calcula igualmente manipulando la función de Lyapunov anterior de manera que la derivada de esta función sea menor que 0. Y al final queda que es equivalente a: gst =1 P1Tz donde >0es un parámetro de diseño y Pes una matriz definida positiva yqueessolucióndeladesigualdadmatricial: P+TP0 8.3 Ley de control: suma de la acción por modelo interno uim y acción estabilizadora ust La acción de control que se implementa es suma de la acción de control por modelo interno ya acción estabilizadora: ub=uim +ust 36
⌧i:SO(3) ⇥R3⇥P!R5 igual al sistema que se ha definido: ⇠i(t)=⌧i(Rod(t),w d od(t),µ(t)) yquecumpleconlaecuacióndemodelointerno,definidaarribayque garantiza la existencia de la acción de control c(w)que hace que nuestro error sea nulo: d dt⌧(Rod(t),w d od(t),µ(t)) = ⌧(Rod(t),w d od(t),µ(t)) c(Rod(t),w d od(t),µ(t)) = ⌧(Rod(t),w d od(t),µ(t)) Por tanto se ha llegado a la conclusión que implementado este modelo para cada uno de los ejes se puede replicar el comportamiento de la acción de control c(w). Replicando lo anterior para los 3 ejes i=1,2,3: ˙ ⇠=⇠ u=⇠ en donde, =0 @ S00 0S0 00S1 A =0 @ Q00 0Q0 00Q1 A Esto se ha implementado en Simulink así: Figure 14: Control modelo interno 43
Como se ve, se ha implementado una integración de la dinámica del modelo interno “ksi”. Y donde el bloque “Dinámica del Modelo Interno Ksi” se ha implementado con el siguiente código: function dksi = fcn(ksi , gst) W= 6 . 2 ∗10^3; %||w(0)|| modulo de la velocidad deseada S=[0 1 0 0 0; 00100; 00010; 00001; 04∗W^ 4 0 5∗W^ 2 0 ] ; phi=[S zeros (5 ,5) zeros (5 ,5); zeros(5,5) S zeros (5 ,5); zeros(5,5) zeros(5,5) S]; dksi = phi∗ksi + gst ; 9.1.3 Estabilidad de la interconexión Después de crear el modelo interno ⇠que contrarresta el efecto del sistema exógeno y añadirlo como una variable más de estado se debe ser asegurar la estabilidad de la conexión (x, ⇠).Paraestoseañadeeltérminogst al modelo interno y que se diseña con argumentos de Lyapunov. El término gst se ha implementado de la siguiente manera : gst =1 P1T(ew+k1eq) function gst = fcn(k1, qerr , werr) beta= 5∗10^3; Q=[1 0 0 0 0 ] ; gamma=[Q z e r o s ( 1 , 5 ) z e r o s ( 1 , 5 ) ; zeros(1,5) Q zeros (1 ,5); zeros(1,5) zeros(1,5) Q]; P= [ 7 . 9 5 , 96.09,694.77,2150.92,4867.73,0,0,0,0,0,0,0,0,0,0; 96.09,1817.74,17243.04,74694.01,112765.08,0,0,0,0,0,0,0,0,0,0; 694.77,17243.04,204788.29,1199232.40,104586.91,0,0,0,0,0,0,0,0,0,0; 2150.92,74694.01,1199232.40,10430890.31,27944049.88,0,0,0,0,0,0,0,0,0,0; 4867.73,112765.08,104586.91,27944049.88,536545695.41,0,0,0,0,0,0,0,0,0,0; 0,0,0,0,0,7.95,96.09,694.77,2150.92,4867.73,0,0,0,0,0; 44
0,0,0,0,0,96.09,1817.74,17243.04,74694.01,112765.08,0,0,0,0,0; 0,0,0,0,0,694.77,17243.04,204788.29,1199232.40,104586.91,0,0,0,0,0; 0,0,0,0,0,2150.92,74694.01,1199232.40,10430890.31,27944049.88,0,0,0,0,0; 0,0,0,0,0,4867.73,112765.08,104586.91,27944049.88,536545695.42,0,0,0,0,0; 0,0,0,0,0,0,0,0,0,0,7.95,96.09,694.77,2150.92,4867.73; 0,0,0,0,0,0,0,0,0,0,96.09,1817.74,17243.04,74694.01,112765.08; 0,0,0,0,0,0,0,0,0,0,694.77,17243.04,204788.29,1199232.40,104586.91; 0,0,0,0,0,0,0,0,0,0,2150.92,74694.01,1199232.40,10430890.31,27944049.88; 0,0,0,0,0,0,0,0,0,0,4867.73,112765.08,104586.91,27944049.88,536545695.41] gst= (1/beta )∗inv(P)∗gamma’ ∗(werr+k1∗[qerr(2);qerr(3);qerr(4)]); Donde P es una matriz de dimensión 15 la cual es solución de la desigualdad matricial que se ha resuelto utilizando el código que viene en el apéndice II. 9.1.4 Simulación del algoritmo de control en Simulink La dinámica del error de la orientación y velocidad angular, el modelo interno ylaaccióndecontroldiseñada,sumadelaacciónpormodelointernoyacción estabilizadora que se simula es: ˙ eq=1 2E(eq)ew J˙ ew=f(R(eq),ew, Rod,w d od,µ)+ub ˙ ⇠=⇠1 P1T(ew+k1eq) ub=⇠k2(1 + new+k1eqn)( ew+k1eq) AcontinuaciónseobservacomosehaimplementadoenSimulink.Seobserva como se ha dicho el bloque de definición de los parámetros, la referencia que actúa como input exógeno en el sistema y a continuación se ve el bucle cerrado de la dinámica del error del satélite y la del controlador. 45
Figure 15: Esquema general del proyecto Acontinuaciónseve,elbloqueenterodelcontroladorcompuestodela acción por modelo interno y la acción estabilizadora: Figure 16: Esquema general del Bloque “Control” Se ve también el bloque de la dinámica del satélite. Como se observa existen dos bloques “tau0” y “taugrav”, estos se corresponden con los momentos de fuerza creados por la inercia del satélite y el efecto de la gravedad. A su vez estos dos bloques se suman y se multiplican por la matriz inversa J para obteniéndose ˙ ewla derivada de la velocidad angular, esta a su vez se integra para hallar ew.Tambiénseobservaladinámicadelaorientacióndel satélite que tiene como input el error de la velocidad angular. A su vez esta se integra para hallar la orientación eq. 46
Figure 17: Esquema general del Bloque “Dinámica del Error del Satélite” 47
Figure 18: Error Orientación del Satélite Los resultados de la simulación se muestran a continuación: 48
Acontinuaciónseobservaelerrordelavelocidadangular: Figure 19: Error de la velocidad angular ew 49
Yseobservalaaccióndecontrol: Figure 20: Acción de control ub(t) 50
9.2 Misión: Vuelo orientado hacia la tierra 9.2.1 Objetivo de la misión Figure 21: Objetivo de la misión: Modo orientado hacia la tierra Análogamente a la misión anterior, las condiciones iniciales del estado de la referencia son: wd od(0) = 0 Rod(0) = I Ylareferenciaaseguirpornuestrosatélitetendrálasiguienteforma: wd od(t)=0 y Rod(t)=I lo cual se trata como una perturbación en el sistema generado por el sistema dinámico: 8 > < > : ˙ Rod =Rod Skew(wd od) ˙wd od =0 ˙µ=0 51
El estado de la incertidumbre en los parámetros µ=✓J w0◆se trata como también como un input exógeno en mi sistema. 9.2.2 Diseño de la acción de control por modelo interno uim Un tipo de misión que se puede simular con este algoritmo de control es el de mantener el satélite en órbita de manera que siempre este observando hacia la tierra. En este caso la orientación del satélite siempre tiene que estar de tal manera que esté en el mismo sistema referencial que el órbital. De esta manera querríamos tener la velocidad angular del satélite nula y la orientación del satélite de tal manera que sea igual al sistema orbital: En este caso las ecuaciones de movimiento del satélite presentadas antes, ya que Rod =Iw d od =0 Y por tanto la velocidad y aceleración angular deseada respecto a la Tierra son: wb ed =e RTwd ed =e RT[wd od w0RT odj]=)wb ed =w0e RTj ˙wb ed =d dt(e RT[wd od w0RT odj]) =)˙wb ed =w0ew⇥(e RTj) La acción de control que se tiene que replicar a través del modelo interno quedaría: ub=c(Rod,w d od,µ)=⌧b 0(I,0,R od,w d od,µ)⌧b grav(I,Rod,µ) Donde el torque inercial es equivalente a: ⌧b 0(I,0,R od,w d od,µ)=w2 0j⇥(Jj)=w2 00 @ Jyz 0 Jxy 1 A Yeltorqueprovocadoporlagravedadserá: ⌧b grav(I,Rod,µ)=3w2 00 @ Jyz Jxz 01 A Ya que todas expresiones son términos constantes. Es decir ni la orientación deseada ni la velocidad angular deseada varía a lo largo del tiempo, se puede concluir que el control de prealimentación que hace que Rob(t)= 52
APÉNDICE II En esta parte se hace el cálculo de la matriz P positiva necesaria para nuestra función de Lyapunov tal que: PA+ATP0 Esto se va a plantear en matlab como la solución P a la desigualdad matricial y que vamos a calcular con el paquete específico de matlab para desigualdas matriciales PA+ATP0 0P<0 El código que se implementa en matlab, es el siguiente: >>% E m p i e z a a q u í , t e n g o c o m o d a t o l a m a t r i z A >> %C o n f i g u r o e l s o l u c i o n a d o r d e LM Is >> s e t l m i s ( [ ] ) >> %A h o r a e s p e c i f i c o l a e s t r u c t u r a d e P >> >> P= l m i v a r ( 1 , [ s i z e ( A , 1 ) 1 ] ) >>% D e f i n o a h o r a e l L M I >> >> l m i t e r m ( [ 1 1 1 P ] , 1 , A , ’ s ’ ) ; >> l m i t e r m ( [ 1 1 2 0 ] , 1 ) ; >> l m i t e r m ( [ 1 2 2 P ] , 1, 1); >> >>% S o l u c i o n o >> >> LMISYS= g e t l m i s ; >> [ t m i n , P s o l ] = f e a s p ( LMISYS ) ; >> P= d e c 2 m a t ( LMISYS , P s o l , P ) 59
References [1] Robust Autonomous GuidanceAlberto Isidori, Lorenzo Marconi and Andrea Serrani [2] Nonlinear Control Systems - Alberto Isidori [3] Nonlinear Control Systems IIAlberto Isidori [4] Introduction to TopologyBert Mendelson [5] Mathematics: Its Content, Methods and MeaningA.D. Aleksandrov, A.N. Kolmogorov, and M.A. Lavrent’ev [6] The Linear Ouptut Regulation ProblemCeSOS-NTNU 2005 - Andrea Serrani [7] Output Regulation for Nonlinear Systems: an OverviewC.I. Byrnes, A. Isidori [8] The early days of geometric nonlinear controlRoger Brockett [9] A remark on the problem of semiglobal nonlinear output regulation [10] Spacecraft Dynamics and Control: An IntroductionAnton H. de Ruiter, Christopher Damaren and James R. Forbes [11] http://www.braeunig.us/space/orbmech.htm [12] https://blog.technavio.com/blog/top-10-satellite-manufacturers-globalspace-industry 60