INPE-12970-TDI/1018 ESTUDO COMPARATIVO DE TÉCNICAS DE CONTROLE DE ATITUDE EM TRÊS EIXOS PARA SATÉLITES ARTIFICIAIS Gilberto Arantes Júnior Dissertação de Mestrado do Curso de Pós-Graduação em Engenharia e Tecnologia Espaciais/Mecânica Espacial e Controle, orientada pelo Dr. Ijar Milagre da Fonseca, aprovada em 23 de fevereiro de 2005. INPE São José dos Campos 2005 629.7.062.2 ARANTES JR, G. Estudo comparativo de técnicas de controle de atitude em três eixos para satélites artificiais / G. Arantes Jr. – São José dos Campos: INPE, 2005. 201p. – (INPE-12970-TDI/1018). 1.Estabilização de veículos espaciais. 2.Controle de atitude. 3.Estabilização em três eixos. 4.rodas de reação. 5.Controle magnético. I.Título. Aprovado (a) pela Banca Examinadora em cumprimento ao requisito exigido para obtenção do Título de Mestrado em Engenharia e Tecnologia Espaciais/Mecânica espacial e Controle Aluno (a): Gilberto Arantes Junior São José dos Campos, 23 de fevereiro de 2005 ”No meio de qualquer dificuldade encontra-se a oportunidade” ALBERT EINSTEIN “A maior recompensa do trabalho do Ser humano não é o que ele(a) ganha com isso, mas sim o que ele(a) se torna com isso.” DESCONHECIDO A meus pais, GILBERTO MOURA ARANTES e CÉLIA DIAS ARANTES, pelo apoio e paciência. E ao tio e amigo MIGUEL BERNARDES DE CASTRO (in memorian). AGRADECIMENTOS A meus pais, pelo amor, compreensão, paciência e por incentivar e acreditar na importância de ir em busca dos sonhos. À minha famı́lia pela certeza de sempre poder contar com seu eterno apoio e incentivo, em especial aos meus Avós, José, Terezinha, Jerônimo e Cacilda. À querida tia Maria José e ao inesquecı́vel tio Miguel (in memorian) por todos os conselhos, incentivos e valiosa torcida. Aos tios Aluisio pela amizade, generosidade e todos os “ensinamentos filosóficos”, e Branca pelo amor de mãe que me foi dedicado, serei eternamente grato. Ao orientador e amigo Prof. Dr. Ijar Milagre da Fonseca pela orientação, diretrizes, conselhos e à valiosa confiança e por acreditar que eu pudesse realizar este trabalho. Ao apoio financeiro dos meus pais e amigos. Ao grande amigo e colega de curso Rolf Vargas por poder dividir minhas dificuldades e pelo seu bom humor, nos momentos difı́ceis. A todos os colegas do curso, pela amizade e companheirismo demonstrados, em especial a Leandro e Cecı́lia, por toda a colaboração. Ao Laboratório LABSIM e todos os técnicos, pela oportunidade de estudos e utilização de suas instalações. Ao Instituto Nacional de Pesquisas Espaciais (INPE), pela oportunidade. Aos professores do INPE em especial a André Fenili, Marcelo Lopes, Waldemar de Castro, Mario Ricci e Evandro Rocco, pelo conhecimento compartilhado. Aos professores Hélio Koiti Kuga e Roberto Vieira da Fonseca Lopes pelas, sugestões na elaboração desse trabalho. A meus amigos da República, mestre Clério, à simpática mexicana Nora Trelles e ao disciplinado Marcos Timbo por toda a paciência e stress compartilhado. À amiga Rose e ao amigo Fred, pelo incentivo. À Coordenação de Aperfeiçoamento de Pessoal de Nı́vel Superior (CAPES), por um ano de bolsa concedida. A todas aquelas pessoas que de uma maneira ou de outra contribuı́ram, e por omissão não constam nessa lista, peço desculpas e agradeço. RESUMO Neste trabalho propõe-se um estudo sobre técnicas de controle de atitude para satélites estabilizados em três eixos e, através da modelagem e simulação computacional, analisar, comparar e fazer um estudo de alternativas/viabilidade de sistemas de controle de atitude (ACS) em três eixos, com base nos requisitos de missões espaciais. Para o estudo de alternativas/viabilidade dos sistemas de controle de atitude foi realizado um estudo comparativo de diferentes técnicas de controle, utilizadas para estabilização de satélites em três eixos, utilizando-se diferentes atuadores, tais como: 1) rodas (de reação e volantes de inércia); 2) bobinas magnéticas. Os procedimentos de estabilização estudados foram: 1) controle de atitude em três eixos utilizando rodas (de reação e volantes de inércia) e bobinas magnéticas; 2) controle de atitude em três eixos utilizando apenas bobinas magnéticas. Foram utilizados a teoria do Regulador Linear Quadrático (LQR) e Regulador Quadrático Gaussiano (LQG) e controladores não lineares baseados em energia (energy based control ) para o desemvolvimento do projeto de controle no modo de estabilização. A teoria do LQR rastreio (tracking) e o controlador Proporcional Derivativo (PD) foram usados no modo de aquisição de atitude. Para a fase de redução da velocidade angular (detumbling) foi utilizado o controlador de Wisniewski ou Bdot. As diferentes configurações dos ACS são discutidas, comparadas e analisadas, visando avaliar o desempenho dos procedimentos de controle aqui desenvolvidos para os modos de detumble, estabilização e aquisição da atitude. A formulação obtida nesse trabalho foi aplicada no controle de atitude do satélite brasileiro EQUARS (Equatorial Atmosphere Research Satellite), que motivou este trabalho. Os resultados obtidos atendem as especificações e os requisitos do satélite. COMPARATIVE STUDY OF THREE AXES ATTITUDE CONTROL TECHNIQUES FOR ARTIFICIAL SATELLITES ABSTRACT This work deals with the study of attitude control techniques for three axes stabilized satellite, and through modeling and computational simulation, analyze, compare and develop a feasibility study for attitude control system based on space missions requirements. The attitude control system feasibility study is carried out by using different control techniques for satellite three-axis stabilization, and different actuators such as 1) reaction wheels and momentum wheels; 2) torque coils. The procedures used for stabilization were: 1) Three axis attitude control by using reaction wheels and momentum wheels combined with torque coils; 2) Three axis attitude control by using torque coils only. Each of these technique was implemented in computer for the attitude control simulations. The control law was based on the Linear Regulator Quadratic (LQR) and Linear Quadratic Gaussian (LQG) for linear systems and on energy based control for no linear systems. These techniques were used for the stabilization mode. The LQR tracking and the proportional derivative (PD) techniques were used for the acquisition mode. The Wisniewski or Bdot approach has been used for detumblig phase. The different configuration results for the control modes are analyzed and discussed in terms of the performance of the control procedures associated with the attitude control modes (detumble, stabilization and acquisition). The control formulation has been applied for the brazilian satellite EQUARS (Equatorial Atmosphere Research Satellite). The results comply with the satellite specification and requirements. SUMÁRIO LISTA DE FIGURAS LISTA DE TABELAS LISTA DE SÍMBOLOS LISTA DE SIGLAS E ABREVIATURAS CAPÍTULO 1– INTRODUÇÃO 31 CAPÍTULO 2– OBJETIVO 37 2.1 Meios e Recursos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.2 Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 CAPÍTULO 3– REVISÃO BIBLIOGRÁFICA 41 3.1 Atitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Controle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3 Literatura do INPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 CAPÍTULO 4– DEFINIÇÕES DE NOTAÇÕES 4.1 51 Sistemas de Referência . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1.1 Referencial Inercial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1.2 Referencial Orbital . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.1.3 Referencial do Satélite . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2 Representação da Atitude . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.2.1 Matriz de Rotação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.2.2 Parâmetros de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2.3 Ângulos de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2.4 Ângulo Eixo Equivalente . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2.5 Rotações Infinitesimais . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3 Sumário de Notações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 CAPÍTULO 5– FORMULAÇÃO DO PROBLEMA 5.1 5.1.1 63 Modelagem Matemática: Cinemática e Dinâmica . . . . . . . . . . . . . . 63 Equações da Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.1.2 Equações da Dinâmica . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.1.3 Torque devido ao Gradiente de Gravidade . . . . . . . . . . . . . . . . 71 5.2 Campo Magnético Terrestre . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.3 Tratamento das Equações da Dinâmica para os Casos Estudados nesse Trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.4 Linearização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.4.1 Linearização do Modelo do Satélite Equipado com Rodas . . . . . . . . 76 5.4.2 Linearização do Modelo do Satélite Equipado com Bobinas . . . . . . . 78 5.5 5.5.1 5.6 Considerações Sobre os Atuadores, os Modelos Matemáticos e Sistemas de Referênica Adotados . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Modelo dos Atuadores . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Torques Ambientais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.6.1 Torque Devido ao Gradiente de Gravidade . . . . . . . . . . . . . . . . 85 5.6.2 Torque Aerodinâmico . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.6.3 Torque de Pressão de Radiação Solar . . . . . . . . . . . . . . . . . . . 88 5.6.4 Torque Devido ao Dipolo Residual . . . . . . . . . . . . . . . . . . . . 88 5.7 Perturbações internas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 CAPÍTULO 6– PROJETO DE CONTROLE 6.1 6.1.1 6.2 91 Modo de Detumble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Controlador Bdot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 Modo de Estabilização . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2.1 Método LQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.2.2 Método LQG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.2.3 Controladores Baseados em Energia . . . . . . . . . . . . . . . . . . . . 101 6.3 Modo de Aquisição . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.3.1 LQR Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.3.2 Proporcional Derivativo . . . . . . . . . . . . . . . . . . . . . . . . . . 104 CAPÍTULO 7– SIMULAÇÕES 109 7.1 Modo de Detumble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 7.2 Modo de Estabilização . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.2.1 Método LQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 7.2.2 Controladores Baseados em Energia . . . . . . . . . . . . . . . . . . . . 124 7.2.3 Método LQG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.3 Modo de Aquisição . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.3.1 LQR Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 7.3.2 Controle PD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 CAPÍTULO 8– CONCLUSÃO 8.1 141 Sugestões para Tabalhos Futuros . . . . . . . . . . . . . . . . . . . . . . 144 REFERÊNCIAS BIBLIOGRÁFICAS 147 APÊNDICE A–CÁLCULO DO GRADIENTE DE GRAVIDADE 157 APÊNDICE B–PROPAGAÇÃO DA ÓRBITA E TRANSFORMAÇÕES 163 B.1 Cálculo da Órbita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 B.1.1 Posicionamento de Satélites - Problema Direto . . . . . . . . . . . . . . 163 B.1.2 Equação de Kepler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 B.1.3 Matriz de Rotação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 B.2 Matriz de Rotação do Sistema Inercial - Sistema do Sátelite . . . . . . . 169 OF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 B.2.1 Matriz RPECI B.2.2 Matriz ROF P OF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 B.2.3 Matriz RBF OF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 B.3 Cálculo da Latitude, Longitude e Altura . . . . . . . . . . . . . . . . . . 175 APÊNDICE C–IMPLEMENTAÇÃO EM SIMULINK 183 APÊNDICE D–PROGRAMAS EM MATLAB 193 D.1 Projeto LQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 D.2 Projeto LQG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 APÊNDICE E– TOOLBOX ATITUDE 199 E.1 Exemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 E.1.1 Equações do Movimento . . . . . . . . . . . . . . . . . . . . . . . . . . 200 E.1.2 Ambiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 LISTA DE FIGURAS 4.1 Sistemas de referência, inercial (ECI), orbital (OF) e do satélite (BF) 4.2 Construção dos ângulos de Euler . . . . . . . . . . . . . . . . . . . . 57 5.1 Seqüência de rotações 3(ψ) − 2(θ) − 1(φ) 5.2 Ilustração do satélite com rodas de reação, bobinas e os referenciais orbital OF (xo , yo , zo ) e do corpo BF (x, y, z) . . . . . . . . . . . . . . 66 5.3 Configuração das bobinas magnéticas no satélite . . . . . . . . . . . . 70 5.4 Campo magnético local Bo usando o modelo IGRF 2000 . . . . . . . 72 6.1 Configuração do controle LQR . . . . . . . . . . . . . . . . . . . . . . 94 6.2 Sistema planta mais controlador . . . . . . . . . . . . . . . . . . . . . 97 6.3 Estrutura básica do sistema de controle LQG 6.4 Estrutura do filtro de Kalman-Bucy . . . . . . . . . . . . . . . . . . . 100 6.5 Configuração do controle LQR tracking . . . . . . . . . . . . . . . . . 104 6.6 Configuração do controlador PD . . . . . . . . . . . . . . . . . . . . . 108 7.1 Velocidade angular do satélite ω bib . . . . . . . . . . . . . . . . . . . . 111 7.2 velocidade angular ω bib . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.3 Dipolo magnético, torque magnético e potência . . . . . . . . . . . . 112 7.4 Dipolo magnético, torque magnético e potência . . . . . . . . . . . . 113 7.5 Ângulos de atitude roll, pitch e yaw em função do tempo para o caso 1116 7.6 Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas para o caso 1 . . . . . . . . . . . . . . . . . . . . . . . . 116 7.7 Velocidades angulares ω bib e ω bob e rotações por minuto das rodas de reação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 7.8 Ângulos de atitude roll, pitch e yaw em função do tempo para o caso 2118 7.9 Torque τ bw , torque de acoplamento e quantidade de movimento angular das rodas para o caso 2 . . . . . . . . . . . . . . . . . . . . . . . . 119 52 . . . . . . . . . . . . . . . 65 . . . . . . . . . . . . . 99 7.10 Velocidades angulares ω bib e ω bob e rotações por minuto das rodas para o caso 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 7.11 Ângulos de atitude roll, pitch e yaw e tempo para o caso 3 . . . . . . 120 7.12 Velocidades angulares ω bib e ω bob e rotações por minuto das rodas para o caso 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 7.13 Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas para o caso 3 . . . . . . . . . . . . . . . . . . . . . . . . 121 7.14 Ganho do controlador em função das órbitas Kc (1, 1) . . . . . . . . . 122 7.15 Ângulos de atitude roll, pitch e yaw em função do número de órbitas, utilizando bobinas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 7.16 Dipolo magnético, torque de controle (τ bm ) e potência das bobinas em função do número de órbitas, utilizando bobinas . . . . . . . . . . . . 124 7.17 Ângulos de atitude em função do número de órbitas . . . . . . . . . . 125 7.18 Dipolo magnético, torque de controle e potência das bobinas em função do número de órbitas . . . . . . . . . . . . . . . . . . . . . . . 126 7.19 Velocidades angulares ω bib e ω bob e campo magnético terrestre local Bb 127 7.20 Ângulos de atitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 7.21 Dipolo magnético, torque de controle e potência das bobinas em função do número de órbitas . . . . . . . . . . . . . . . . . . . . . . . 129 7.22 Velocidades angulares ω bib e ω bob e campo magnético local Bb em função do número de órbitas . . . . . . . . . . . . . . . . . . . . . . . 129 7.23 Ângulos de atitude estimados em função do tempo . . . . . . . . . . . 131 7.24 Ângulos de atitude simulados . . . . . . . . . . . . . . . . . . . . . . 131 7.25 Ângulos de atitude simulados e estimados . . . . . . . . . . . . . . . 132 7.26 Erro dos ângulos de atitude e variação da atitude . . . . . . . . . . . 132 7.27 Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 7.28 Ângulos de atitude roll, pitch e yaw em função do tempo para a aquisição de atitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 7.29 Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas em função do tempo . . . . . . . . . . . . . . . . . . . . 136 7.30 Rotações por minuto das rodas de reação e velocidades angulares ω bib e ω bob . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 7.31 Ângulos de Euler em função do tempo . . . . . . . . . . . . . . . . . 138 7.32 Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas em função do tempo . . . . . . . . . . . . . . . . . . . . 139 7.33 Rotações por minuto das rodas e velocidades angulares ω bib e ω bob em função do tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 B.1 Referenciais (inercial, da órbita) e os elementos Keplerianos (i, ω, Ω) . 164 B.2 Referencial inercial (ECI) e pseudo orbital (POF) . . . . . . . . . . . 170 B.3 Orientação do referencial orbital (OF) e do referencial pseudo orbital (POF) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 B.4 Referencial inercial e cartesiano geocêntrico . . . . . . . . . . . . . . . 179 B.5 Longitude (λ), latitude (φ) e altura (h) do satélite . . . . . . . . . . . 179 C.1 Implemetação em SIMULINK do Modo de detumble . . . . . . . . . 183 C.2 Lei de controle Bdot . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 C.3 Implementação em SIMULINK do modo de estabilização: satélite equipado com bobinas . . . . . . . . . . . . . . . . . . . . . . . . . . 184 C.4 Controle LQR e controladores baseados em energia . . . . . . . . . . 185 C.5 Modo de estabilização: satélite equipado com rodas, usando o método LQR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 C.6 Modo de estabilização: satélite equipado com rodas, usando o método LQG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 C.7 Modo de aquisição usando o controlador PD . . . . . . . . . . . . . . 186 C.8 Lei de controle LQR tracking . . . . . . . . . . . . . . . . . . . . . . 187 C.9 Lei de controle PD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 C.10 Modelo dinâmico do satélite equipado com bobinas . . . . . . . . . . 187 C.11 Modelo dinâmico do satélite equipado com rodas . . . . . . . . . . . . 188 C.12 Modelo das rodas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188 C.13 Modelo das bobinas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189 C.14 Modelo do torque de gradiente de gravidade . . . . . . . . . . . . . . 189 C.15 Modelo do campo magnético . . . . . . . . . . . . . . . . . . . . . . . 190 C.16 Propagação da órbita . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 C.17 Cálculo da latitude, longitude e altura . . . . . . . . . . . . . . . . . 191 C.18 Matriz de transformação entre os sistemas de referência ECI e BF . . 191 E.1 Bloco Gyrostat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200 E.2 Bloco gradiente de gravidade . . . . . . . . . . . . . . . . . . . . . . . 201 LISTA DE TABELAS 5.1 Elementos orbitais do satélite EQUARS . . . . . . . . . . . . . . . . 73 7.1 Parâmetros de simulação . . . . . . . . . . . . . . . . . . . . . . . . . 110 7.2 Condições iniciais da simulação para o modo de estabilização utilizando a metodologia LQR . . . . . . . . . . . . . . . . . . . . . . . . 114 7.3 Parâmetros de simulação . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.4 Condições iniciais da simulação para o modo de aquisição utilizando rodas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 8.1 Alternativas para o ACS . . . . . . . . . . . . . . . . . . . . . . . . . 141 B.1 Programa MATLAB para o cálculo da órbita circular . . . . . . . . . 169 BF B.2 Programa MATLAB para o cálculo de c1i , (i = 1, 2, 3) da matriz RECI 176 BF B.3 Programa MATLAB para o cálculo de c2i , (i = 1, 2, 3) da matriz RECI 177 BF B.4 Programa MATLAB para o cálculo de c3i , (i = 1, 2, 3) da matriz RECI 178 B.5 Programa MATLAB para o cálculo da longitude, latitude do ponto sub-satélite e altura do satélite . . . . . . . . . . . . . . . . . . . . . . 181 D.1 Programa MATLAB para o projeto LQR do modelo satélite equipado com rodas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194 D.2 Programa MATLAB para o projeto LQR do modelo satélite equipado com rodas (cont.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 D.3 Programa MATLAB para o projeto LQR do modelo satélite equipado com bobinas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 D.4 Programa MATLAB para o projeto LQG do modelo satélite equipado com rodas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 D.5 Programa MATLAB para o projeto LQG do modelo satélite equipado com rodas (cont.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 LISTA DE SÍMBOLOS Bb Bo In J Ji Jo Js Jw Kd Kp L Lb Qc Qf Rab Rc Rf Sa T A B C G(s) Ho Jx , Jy , Jz K(s) Kc Kf Pc Pf To V X, Y, Z x̂ cb1 cb2 cb3 fw hs Vetor campo geomagnético expresso no sistema BF Vetor campo geomagnético expresso no sistema OF Matriz identidade de ordem n Tensor de inercial em relação ao sistema BF Índice de desempenho Tensor de inercial em relação ao sistema o Tensor de inercial do satélite sem as rodas Tensor de inercial das rodas Matriz de ganhos derivativos do controlador PD Matriz de ganhos proporcionais do controlador PD Quantidade de movimento angular do satélite referido e expresso no sistema inercial Quantidade de movimento angular do satélite referido ao sistema inercial e expresso no sistema BF Matriz de ponderação de estados Matriz de covariânça do ruı́do nas medidas Matriz de rotação do sistema b para o sistema a Matriz de ponderação da lei de controle Matriz de covariânça do ruı́do da dinâmica Operador anti-simétrico Torque interno Matriz do sistema Matriz dos atuadores Matriz de sensores Matriz função de transferência de um sistema Quantidade de movimento angular nominal no eixo de pitch Momentos principais de inércia Função de transferência do controlador Matriz de ganhos do regulador LQR Matriz de ganhos do filtro de Kalman Matriz solução da equação de Riccati no regime estacionário para o caso do regulador Matriz de covariânças dos estados estimados Perı́odo orbital Função de Lyapunov do sistema Coordenadas do sistema inercial Estimativa do estado x Cosenos diretores de xo em relação aos eixos do corpo x, y, z Cosenos diretores de yo em relação aos eixos do corpo x, y, z Cosenos diretores de zo em relação aos eixos do corpo x, y, z Sinal de comando Quantida de movimento angular do satélite hw lw mb mc q r A D e e n gm , hnm h h i ik k N u v w x, y, z xo , yo , zo y Ω Φ β ω bbw ω bib ω biw ω bob τ bbc τ bcw τ bc τ bd τ bg τ bm τa τ ext τs Quantida de movimento angular das rodas relativo ao satélite Quantida de movimento angular das rodas Momento dipolo magnéico das bobinas Dipolo magnético residual para o detumble Quaternion Sinal de referência Área da secção da bobina magnética Matriz de influênça do controle na saı́da Excentricidade da órbita Sinal de erro Coeficiente de Gauss Altura do satélite Ganho do controlador EBC com realimentação de velocidade angular Inclinação da órbita Corrente que passa pela bobina magnética disposta na direção k Ganho do controlador Bdot Número de espiras Sinal de contole Ruı́do branco na medida Ruı́do branco na dinâmica Coordenadas do sistema do satélite Coordenadas do sistema orbital Sinal de saı́da do sistema Ascensão do nodo ascendente Ângulo de rotação em torno do vetor unitário λ Ganho do controlador EBC com realimentação de atitude Vetor do quaternion Vetor velocidade angular das rodas referido ao sistema BF, expresso no sistema BF Vetor velocidade angular do satélite referido ao sistema ECI, expresso no sistema BF Vetor velocidade angular das rodas referido ao sistema ECI, expresso no sistema BF Vetor velocidade angular do satélite referido ao sistema OF, expresso no sistema BF Torques de controle das bobinas magnéticas ou magneto torques Torques de controle das rodas Torques de controle expresso no sistema BF Torques extenos que agem sobre o satélite expresso no sistema BF Toque de gradiente de gravidade Torque magnético devido ao dipolo residual Toque de arraste aerodinâmico Torque(s) exteno(s) que agem sobre o satélite Toque de pressão de radiação solar η ωo φr , θr , ψr φ ψ θ Parte real do quaternion Velocidade orbital média Atitude de referência Ângulo de rotação em torno do eixo de roll Ângulo de rotação em torno do eixo de yaw Ângulo de rotação em torno do eixo de pitch LISTA DE ABREVIATURAS E SIGLAS ACS BF EBC ECI EQUARS FG IAGA IGRF INPE LQG LQR MIMO OF OOF PD PMM POF SISO VLHL Sistema de Controle de Atitude Referencial do Satélite Controladores Baseados em Energia Referencial Inercial Equatorial Atmosphere Research Satellite Rereferencial Cartesiano Terrestre Association of Geomagnetism and Aeronomy International Geomagnetic Reference Field Istituto Nacional de Pesquisas Espaciais Regulador Quadrático Gaussiano Regulador Linear Quadrático Multiplas Entradas Multiplas Saı́das Referencial Orbital Referencial da Órbita Proporcional Derivativo Plataforma Multi-Missão Referencial Pseudo Orbital Única Entrada Única Saı́da Vertical Local e Horizontal Local CAPÍTULO 1 INTRODUÇÃO A orientação de um satélite, em relação a um sistema de referência conhecido, é denominada atitude e o movimento de rotação em torno do seu centro de massa é denominado movimento de atitude. Assim, a atitude e o movimento de atitude especificam a orientação espacial e o movimento rotacional em torno do centro de massa do satélite (Wertz, 1978 e Hughes, 1986). Para determinar a atitude de um satélite em relação a um sistema de referência, esse deve estar equipado com sensores que possam fornecer a sua orientação em relação ao Sol, a Terra, e/ou alguma estrela fixa bem como em relação ao vetor campo magnético terrestre (direção e magnitude). Em certos casos é necessário considerar elementos orbitais do satélite para que seja possı́vel determinar completamente a atitude e o movimento de atitude do veı́culo espacial. A análise de atitude pode ser dividida em determinação, previsão e controle de atitude (Wertz, 1978). Este trabalho trata da previsão e o controle de atitude. A previsão é o processo de prever a orientação do satélite pelo uso de modelos que permitam extrapolar sua atitude, sendo necessário conhecer as forças perturbadoras que agem sobre o satélite, portanto, necessário modelá-las. O controle é o processo de orientar o satélite de maneira que esse adquira ou mantenha a atitude prefixada pela missão. A implementação do modelo cinemático e dinâmico do satélite é feita neste trabalho utilizando o software MATLAB e SIMULINK, bem como a implementação das leis e estratégias de controle de atitude em três eixos. Satélites artificiais, em uma grande variedade de missões espaciais, seja para fins meteorológicos, de telecomunicações, de sensoriamento remoto e cientı́fico, devem estar com uma ou mais faces/equipamentos orientada(o)s para direções especı́ficas, tais como para a Terra (BrasilSAT, Syncom, etc), para o Sol (Soho) ou para outras estrelas (Hubble). É evidente a necessidade de se conhecer a atitude do satélite, de modo que se possa estabilizar seu movimento de acordo com a atitude nominal especificada. Isto é feito através de um sistema de controle, projetado de acordo com os requisitos da missão. O controle de atitude é, portanto, o processo de orientar o satélite de maneira que este adquira ou mantenha a atitude nominal especificada. Este processo pode ser visto sob dois modos: modo de aquisição de atitude, que 31 consiste em levar o satélite para a atitude nominal a partir de uma atitude qualquer e modo normal de operação, que consiste em estabilizar/manter a atitude, de acordo com a atitude nominal, fazendo correções, quando necessárias. Durante a fase de aquisição de atitude são, em geral, necessárias manobras de grandes ângulos, uma vez que o sistema de controle de atitude (ACS) deve ser capaz de, a partir de uma atitude qualquer, levar o satélite para a atitude nominal. Já a fase de estabilização/manutenção, em geral, requer manobras de pequenos ângulos objetivando correções na atitude. Eventualmente os requisitos do ACS podem impor manobras de grandes ângulos para reorientar o satélite em determinada direção. Um requisito deste tipo requer manobra de reaquisição da atitude. As duas configurações básicas para estabilizar satélites são aquelas em um eixo (exemplos do SCD1 e SCD2) e em dois eixos, essa conhecida como estabilização em três eixos (exemplos do CBERS1 e CBERS2), uma vez que estabilizando-se dois eixos garante-se a estabilização do terceiro eixo. O estudo da estabilização em um eixo é descrito em Kuga L. D. e Guedes (1987), Kuga e Guedes (1987), Quirelli (2002) e Zanardi et al. (2003a). Esse trabalho é dedicado ao estudo da estabilização em três eixos de satélites artificiais. A estabilização em três eixos pode ser descrita em relação ao sistema orbital. Nesse referencial o movimento em torno da direção da velocidade orbital é denominado roll (rolamento). O movimento em torno da direção normal à órbita é denominado pitch (arfagem), e finalmente o movimento em torno da direção Nadir/Zênite é denominado yaw (guinada). Satélites estabilizados em três eixos podem ser vistos como um sistema com quantidade de movimento angular nulo (zero momentum system) ou então como um sistema em que a quantidade de movimento angular é não nula (bias momentum system). Os sistemas com quantidade de movimento angular não nulo contém, geralmente, um volante de inércia com quantidade de movimento angular nominal diferente de zero (Wertz, 1978 e Wie, 1998). Este equipamento difere das rodas de reação pelo fato de ter velocidade angular nominal diferente de zero em certa direção (Larson e Wertz, 1992). Atualmente o programa espacial brasileiro desenvolve, entre outros, projetos de dois satélites; o Satélite Equatorial Upper Atmosphere Research Satellite de aplicações cientı́ficas (EQUARS) e a plataforma multi-missão (PMM), visando uma plataforma única para satélites com diferentes missões. Os requisitos dessas missões impõem um sistema de controle em três eixos. Essas duas missões espaciais consti32 tuem a principal motivação para o desenvolvimento do trabalho. Diversas técnicas de estabilização podem ser utilizadas quando se deseja estabilizar satélites em três eixos. Dentre elas destacam-se nesse trabalho, aquelas que empregam como atuadores: rodas de momentum e bobinas magnéticas. Nesse trabalho são estudados/analisados três modos de operação: • Detumble: reduzir a velocidade angular do satélite. • Estabilização: adquirir e manter o referencial do satélite alinhado com o referencial orbital (vertical local e horizontal local). • Aquisição: aquisição de uma orientação arbitrária. Os modelos matemáticos do veı́culo espacial usados nos procedimentos de controle de atitude avaliados para o modo de detumble e o modo de estabilização são: rodas de reação combinadas com bobinas magnéticas e um sistema empregando somente bobinas magnéticas. As bobinas são utilizadas para a redução da velocidade angular, ou seja, no modo de detumble e desaturação. Entretanto, nesse estudo a estratégia de dessaturação não é estudada. Para o modo de aquisição é avaliado o emprego de rodas com objetivo de analizar a exequibilidade do procedimento a partir das especificações das rodas para o satélite EQUARS. As leis e estratégias de controle empregadas nesse trabalho são associadas à técnica do Regulador Linear Quadrático (LQR), LQR tracking, Regulador Quadrático Gaussiano (LQG), controle proporcional derivativo (PD), controladores baseados em energia (attitude feedback e angular velocity feedback ) usando a teoria de Lyapunov e o controlador Bdot que também é baseado na teoria de Lyapunov. As estratégias de controle adotadas para os três modos de operação são: • O controlador Bdot, usando apenas medidas dos magnetômetros, para o descapotamento/detumbling do satélite. Os atuadores empregados são as bobinas magnéticas. • Os controladores LQR/LQG e os baseados em energia para o modo de estabilização. A metodologia LQR é usada para rodas e bobinas e a metodologia 33 LQG apenas para rodas. Os controladores baseados em energia são usados com o emprego de bobinas magnéticas. • A metodologia LQR tracking e o controlador PD para o modo de aquisição, com o objetivo de avaliar a exequibilidade/factibilidade do emprego das rodas especificadas para o satélite brasileiro EQUARS. A dissertação está organizada de acordo com a seqüência descrita a seguir. No Capı́tulo 1 é feita uma breve introdução sobre o trabalho, descrevendo brevemente os procedimentos e estratégias de controle adotadas. Nesse Capı́tulo é incluı́da a motivação do estudo, que se refere à aplicação desse estudo ao satélite brasileiro EQUARS. O Capı́tulo 2 apresenta os objetivos do trabalho. O Capı́tulo 3 apresenta a revisão bibliográfica sobre o tema. No Capı́tulo 4 é realizada uma revisão das definições e fundamentos matemáticos envolvendo atitude e movimento de atitude de veı́culos espaciais. No Capı́tulo 5 são apresentados os modelos matemáticos da cinemática e dinâmica do veı́culo espacial. É mostrado o modelo do campo magnético terrestre e o modelo dos principais torques ambientais (gradiente de gravidade, pressão de radiação, arrasto atmosférico e dipolo residual). No Capı́tulo são discutidos brevemente os torques de origem interna e uma aproximação do seu modelo. É apresentado ainda a linearização dos modelos dinâmicos e cinemáticos da atitude para os diferentes procedimentos adotados. No Capı́tulo 6 são apresentados as estratégias de controle, os fundamentos e metodologias usadas na formulação das leis de controle empregadas nos sistemas de controle de atitude (ACS). No Capı́tulo 7 são apresentadas as simulações da dinâmica e controle do veı́culo para as diferentes configurações do ACS, envolvendo os diferentes modos de operação. O estudo é aplicado ao satélite brasileiro EQUARS. É feita a análise dos resultados 34 e comparações das estratégias de controle empregadas. Capı́tulo 8 encerra-se com a análise e discussões dos resultados, sendo feitas sugestões para trabalhos futuros. 35 36 CAPÍTULO 2 OBJETIVO Este trabalho visa modelar, simular, analisar, comparar resultados e apresentar um estudo de alternativas de sistemas de controle de atitude (ACS) para satélites estabilizados em três eixos em diferentes modos de operação, utilizando como atuadores: • Rodas de reação/volantes de inércia (reaction wheel /momentum wheel ); • Bobinas magnéticas. Os procedimentos adotados para os modos de detumble e estabilização usando os atuadores citados acima são: 1) satélite equipado com rodas e bobinas magnéticas; 2) satélite equipado apenas com bobinas magnéticas. O projeto de controle usado para o modo de estabilição utiliza diferentes reguladores: 1) os baseados nas técnicas modernas de controle multivariável: Regulador Linear Quadrático (LQR) e Regulador Quadrático Gaussiano (LQG); 2) os controladores baseados em energia (energy based control ) com realimentação de velocidade angular (angular velocity feedback ) e realimentação de atitude (attitude feedback ), ambos utilizando a teoria de Lyapunov. O projeto de controle usado para o modo de detumble usa apenas medidas dos magnetômetros e é baseado na teoria de Lyapunov. O projeto de controle para o modo de aquisição utiliza: 1) controlador Proporcional Derivativo (PD); 2) LQR rastreio/tracking. Esse projeto tem o objetivo de avaliar e discutir a factibilidade do emprego das rodas especificadas para o satélite brasileiro EQUARS, em adquirir uma atitude arbitrária, e ainda comparar os controladores usados. O objetivo principal desse trabalho é realizar um estudo de alternativas e um estudo comparativo em função do desempenho e análise de custo & benefı́cio dos diferentes sistemas de controle de atitude. As configurações dos ACS(s) estudados serão aplicados ao satélite brasileiro EQUARS (Equatorial Upper Atmosphere Research Satellite). O estudo desenvolvido nesse trabalho (estudo de alternativas/viabilidade) pode se extendido por exemplo, a plataforma Multi-Missão (PMM). 37 2.1 Meios e Recursos Os principais meios e recursos utilizados para se alcançar os objetivos propostos no trabalho são: • pesquisa bibliográfica utilizando os recursos da biblioteca do INPE, artigos publicados na literatura da área, e pesquisa no portal de periódicos da CAPES; • uso dos recursos da divisão de Mecânica Espacial e Controle, no que se refere a meios computacionais para simulação, laboratório LABSIM, incluindo os ambientes MATLAB e SIMULINK para o desenvolvimento das simulações. 2.2 Metodologia O desenvolvimento do trabalho proposto está fundamentado nas seguintes metodologias: • Modelagem matemática da dinâmica de atitude para um veı́culo espacial contendo como atuadores: rodas (de reação e/ou volantes de inércia) e bobinas magnéticas, incluindo os efeitos do gradiente de gravidade. • Análise dinâmica do modelo matemático, incluindo a linearização do modelo, usados no desenvolvimento do projeto de controle (reguladores), para os controladores lineares (PD, LQR e LQG). • Modelagem e implementação, no ambiente SIMULINK, da órbita e do campo magnético terrestre, modelo IGRF (International Geomagnetic Reference Field ). • Estudo do problema do regulador linear quadrático (LQR) e regulador quadrático gaussiano (LQG). • Estudo dos controladores baseados em energia, utilizando a teoria de Lyapunov. • Estudo do controlador proposto por Wisniewski (1996) ou Bdot. 38 • Implementação dos sistemas dinâmicos no ambiente MATLAB e SIMULINK. • Simulação dos sistemas dinâmicos (não lineares). • Análise e comparação das estratégias de controle empregadas nas diferentes configurações dos ACS(s), desenvolvidos para satélites estabilizados em três eixos, nos diversos modos de operação. Aplicação do estudo ao satélite brasileiro EQUARS. 39 40 CAPÍTULO 3 REVISÃO BIBLIOGRÁFICA Este Capı́tulo apresenta uma revisão da literatura relacionada ao principais assuntos envolvidos nesse trabalho: controle de atitude de satélites, técnicas de estabilização em três eixos, modelamento matemático dos sistemas (satélite, atuadores, meio) e estratégias de controle. As metodologias modernas de controle multivariável LQR, LQR/tracking e LQG são revistas, controladores PD e reguladores não lineares baseados em energia também são apresentados. 3.1 Atitude As referências Wertz (1978), Kaplan (1976), Wie (1998) e Hughes (1986) apresentam um estudo do movimento de atitude de satélites artı́ficiais, mostrando o histórico, os fundamentos e os conceitos fı́sicos fundamentais utilizados na previsão, controle e determinação de atitude. A modelagem, análise e controle de sistemas dinâmicos são revistos, as técnicas de controle de atitude para satélites artificiais estabilizados em três eixos também são discutidas. Os tipos e descrição dos sensores e atuadores empregados no movimento de atitude são apresentados em Wertz (1978) e em Pilchowski (2001). Uma descrição do uso de sensores em satélites é mostrada em Wright e Wong (1989). As várias formas de representação/parametrizações da atitude (quaternions, ângulo eixo equivalente, ângulos de Euler, cossenos diretores, variáveis de Andoyer, parâmetros de Gibbs) são apresentadas em Wertz (1978), Rodrigues e Zanardi (2004), Wie (1998), Fauske (2003) e Junkins e Turner (1986). O problema de estabilização em três eixos é também discutido em Martins Neto (2001). A análise de missão é uma fase na qual define-se o que deve ser feito sem necessariamente definir como fazer. A análise de missão envolvendo especificações do sistema de controle de atitude para estabilização em três eixos é descrita e discutida em Larson e Wertz (1992). Em Marteau e Rogers (1996) é feita uma análise de sistemas de controle de atitude (ACS), abordando aspectos gerais sobre as especificações (atuadores, sensores, computadores de bordo) de um sistema de controle de atitude (ACS). O Autor discute um ACS com custo “razoável”, para pequenos satélites (< 500kg) que requeiram apontamento de “razoável”precisão (20 arcsec), precisão tı́pica requerida para satélites já operacionais como IRAS (Infrared Astronomical 41 Satellite) (Beichman et al., 2004), ASCA (Advanced Satellite for Cosmology and Astrophysics) (Tamura, 1998) e SOHO (Solar & Heliospheric Observatory) (Gurman, 2004). A modelagem da dinâmica de atitude utilizando-se rodas (de reação e volantes de inércia) para estabilização em três eixos é apresentada nas referências Wie (1998) (onde se utilizam duas rodas de reação e um volante de inércia), Yairi (1994) e Fichter e Zentcraf (1996). Yairi (1994) e Fichter e Zentcraf (1996) apresentam a mesma modelagem de Wie (1998), mas utilizam quatro rodas de reação, sendo uma delas disposta na diagonal (skew symmetric), redundante. A roda redundante é disposta de tal forma que forneça uma quantidade de movimento angular nas três direções principais de inércia, o que pode, eventualmente, em caso de falhas substituir uma ou mais rodas ao longo dessas direções. Essas referências mostram o desenvolvimento das equações do movimento rotacional de um corpo rı́gido equipado com rodas. O problema ótimo no procedimento de estabilização de satélites artifı́ciais com a utilização de rodas é apresentado em El-Gohary (2003). Este autor apresenta um estudo de estabilidade segundo Lyapunov. Em Varatharajoo e Fasoulas (1975) é feita a análise do problema de atitude empregando rotores para satélites de observação da Terra. Spindler (2000) aborda o problema de controle para estabilização em três eixos empregando N volantes de inércia, exemplificando o procedimento para o satélite MONS-ballerina. O problema de estabilização de atitude em três eixos usando apenas atuadores eletromagnéticos é discutida em Kaplan (1976), Wertz (1978), Bushenkov e Smirnov (2002), Psiaki (2001), Wang e Shtessel (1998), Wisniewski e Blanke (1999) e Musser e Ebert (1989). Em Wertz (1978) e Carrara (1982) é encontrada a modelagem das perturbações ambientais (forças e torques) que atuam sobre o satélite no espaço, devido ao campo magnético, o campo gravitacional e a radiação proveniente do Sol e da Terra. A modelagem do campo magnético, modelo IGRF (International Geomagnetic Reference Field ) usado nesse trabalho é encontrado em Macmillan e Quinn (2000). Outros modelos do campo magnético (dipolo e quadripolo) podem ser encontrados em Zanardi et al. (2003b) e Zanardi et al. (2004). O Desenvolvimento das equações da dinâmica de atitude usando atuadores eletromagnéticos/bobinas magnéticas é apresentado em Psiaki (2001), Musser e Ebert (1989), Wisniewski (1997) e Marteau e Psiaki (1988). Outras referências também 42 trazem as equações da dinâmica e cinemática de atitude, como Bushenkov e Smirnov (2002), Fauske (2002), Grassi e Moccia (1995) e Wang e Shtessel (1998). Essas referências discutem a utilização de bobinas magnéticas para o controle para pequenos satélites visando menor custo para precisão requerida de 1o a 2o , sem requerer alta fonte de energia. Wisniewski e Blanke (1999) analisam a técnica para o satélite dinamarquês Orsted (Clausen, 2004) (satélite aplicado no estudo do campo magnético da Terra). O autor mostra que é factı́vel obter estabilização em três eixos utilizando apenas torques magnéticos para satélites de órbita baixa (LEO), sujeitos ao gradiente de gravidade. O ASRI (Australian Space Research Institute) (Ardebil, 2004) confirma a viabilidade de se usar apenas atuadores magnéticos para o controle de atitude de pequenos satélites, através da missão cientı́fica TechSAT. Silani e Lovera (2003) apresentam uma revisão do problema de estabilização de atitude em três eixos para pequenos satélites, usando atuadores magnéticos (bobinas), baseado na teoria de controle linear e não linear. Junkins e Carrington (1980) realizam um estudo de otimização para manobras de atitude usando atuadores magnéticos. A referência Kim e Choi (1999) trata da utilização de um sistema de controle de atitude (ACS), utilizando três rodas de reação combinadas com bobinas magnéticas para o microsatélite Koreano KITSAT-3 (satélite de telecomunicações), visando estabilização em três eixos, em órbitas baixas (LEO). Nos requisitos de payload /carga útil do KITSAT-3 são exigidas alta precisão de apontamento (0.05o ) e estabilidade (0.014o rad/s) da plataforma. Os atuadores magnéticos também podem ser utilizados para a desaturação das rodas. Bang e Choi (2003) analisam a desaturação de rodas utilizadas para manobras de grandes ângulos. Outra técnica de controle para estabilização em três eixos, desenvolvida primeiramente para os satélites TIROS (Television Infrared Observation Satellite), por Harold Perkel, utiliza uma configuração particular de um volante de inércia (momentum wheel ), ao longo do eixo de arfagem (pitch), combinando com dois atuadores eletromagnéticos, dispostos ao longo do eixo de rolamento (roll ) e guinada (yaw ). Esta técnica é denominada stabilite e é encontrada na referencia Perkel (1966). Em Hamzah e Hashida (1999) é descrito o desenvolvimento do sistema de controle de atitude para o satélite TiungSAT-1, utilizando-se a configuração de Perkel obtendo uma precisão de ampontamento de ±1o . Whitford e Forrest (1998) descrevem o desenvolvimento de um sistema de controle de atitude (ACS) para estabilização em três eixos também baseado na combinação 43 de rodas (de reação e volantes de inércia) e bobinas magnéticas. Os vários modos de operação são avaliados para sistema de controle de atitude (ACS). O ACS é aplicado ao satélite CATSAT (Co-operative Astrophysical and Technology Satellite) do programa STEDI (Student Explorer Demonstration Initiative). Um dos principais problemas inerentes à utilização de rodas é o da saturação no processo de controle, devido a manobras de grandes ângulos (Bang e Choi, 2003) e/ou torques seculares. Em Buckingham e Smirnov (1972) é discutida a desaturação de rodas de reação empregando atuadores magnéticos. Em Gökçev e Meerkov (2001) é apresentada a metodologia para sistemas com saturação de atuadores utilizados na estabilização de satélites sujeito a perturbações seculares. Da vasta literatura que apresenta o desenvolvimento das equações da dinâmica não linear, para satélites equipados com bobinas, destacam-se os trabalhos de Fauske (2002), Cohen (1973), Spencer (1977), Shigehara (1972) e Alfriend (). As equações do sistema tornam-se variantes no tempo devido ao campo geomagnético. A análise da controlabilidade de satélites equipados apenas com bobinas magnéticas é encontrada em Bhat e Dham (2003). Wisniewski e Markley (1999) desenvolvem uma metologia de controle ótimo aplicada ao controle de atitude em três eixos. Propriedades elétricas de materiais utilizados em bobinas magneticas são apresentados e discutidos em Legg (2003). A descrição, especificações e precisões de rodas (de reação e volantes de inércia) da TELDIX utilizadas para estabilização em três eixos são encontrados em Auer (1983) e Heidelberg (2004). Em Auer (1983) são também apresentadas repostas para os comandos de controle. 3.2 Controle Nesta seção é feita um revisão bibliográfica das metodologias LQR, LQR rastreio (tracking) e LQG. Uma revisão de controladores baseados em energia, também reguladores, obtidos a partir da teoria de Lyapunov é apresentada. O controlador usado para o modo de detumble, conhecido como Bdot é revisto. Metodologia LQR e LQG A teoria do regulador linear quadrático (LQR) e do regulador linear gaussiano (LQG) aplicadas em controle multivariável é encontrada em uma vasta literatura. As prin44 cipais referências estudadas pelo autor para elaboração e implementação dessas metodologias no sistema de controle de atitude foram Dorato e Cerone (1995), Maciejowski (1989), Kwakernaak e Sivan (1972), Moore e Anderson (1990), Brown (1997) e Kirk (1970). Essas referências fornecem os métodos para formular a lei de controle para sistemas dinâmicos. O projeto de controle LQR e LQG é baseado na linearização dos sistemas dinâmicos, definindo uma função objetivo a ser minimizada, e na obtenção de uma matriz de ganhos (variantes no tempo ou não) usada na realimentação (Overby, 2004), o método se aplica a sistemas variantes no tempo, como o caso de satélites equipados com bobinas magnéticas, tornando-se uma alternativa atraente, devido a robustez e fácil implementação, essa metodologia é usada por Wisniewski e Markley (1999), Wisniewski (1996) e Marteau e Psiaki (1988) para o sistema de controle de atitude de satélites equipados com atuadores magnéticos. Yairi (1994) usa a formulação LQR para o sistema de controle de atitude de satélites estabilizados em três eixos equipados com rodas. Nesse trabalho a lei de controle, formulada a partir do sistema linear, é aplicada ao sistema dinâmico completo, não linear, assim como feito por Musser e Ebert (1989) e Fauske (2002). No controle de atitude para satélites estabilizados em três eixos é necessário manter a atitude em relação a um sistema de referência, inercial ou não inercial, como é o caso de satélites Terra apontados em que o sistema do satélite deve estar alinhado com a vertical local e horizontal local (VLHL) ou referencial orbital, um sistema de referência móvel (Wertz, 1978). Portanto, no problema de estabilização a referência para o controlador é zero ou constante, tratando-se de um regulador (Distefano e Willians, 1972). Em Kirk (1970) e Dorato e Cerone (1995) é apresentada a formulação das lei de controle para o caso em que a referência difere de zero, essa formulação é conhecida como LQR rastreio/tracking (Dorato e Cerone, 1995). Em Overby (2004) é discutido várias estratégias de controle e reguladores usados para a estabilização do satélite Norueguês Ncube. Entre essas estratégias a metodologia LQR é apresentada, mostrando boas propridades de estabilidade. Overby (2004) e Kristiansen (2000) apresentam um algoritmo para a seleção dos ganhos das matrizes de ajuste usadas no LQR, baseado na saturação dos atuadores e intervalos de operação para os estados. Em Dorato e Cerone (1995), Kwakernaak e Sivan (1972) e Souza (1987) é discutida a escolha/seleção das matrizes dos ganhos com base nos requisistos e compromissos/trade-off do processo. Na metologia LQR não 45 é levada em conta as incertezas do problema, ou seja, a metodologia LQR é determinı́stica. Através da utilização da metodologia LQG, um controlador baseado em observador, essas incertezas (estocásticas) podem ser atribuı́das à planta ou processo. A solução do problema do LQG é obtida pelo uso do princı́pio da separação que possibilita a separação do problema original (estocástico) em dois problemas: 1) estimação ótima do estado, de modo que sua covariância seja minimizada. Esse sub-problema é resolvido pelo uso de um filtro de Kalman, ignorando-se completamente o problema de controle; 2) formulação da lei de controle usando a metodologia LQR, fazendo uso da estimativa como se fosse a medida exata do estado, ignorando-se completamente os aspectos ou natureza estocásticas do processo (Flora, , Moore e Anderson, 1990 e Dorato e Cerone, 1995). A solução do problema 1 consiste na determinação da matriz ganho ou ganho de Kalman, a estrutura do filtro é encontrada em Maciejowski (1989), Kwakernaak e Sivan (1972) e Athans (1971). A solução do problema 2 consiste na derminação da matriz ganho do controlador usado na realimentação. Para a solução do problema LQG ser assintoticamente estável é necessário que o sistema seja completamente observável, para a solução do problema 1 e completamente controlável, para a solução do problema 2 (Brown, 1997 e Dorato e Cerone, 1995). Controladores baseados em energia Uma importante ferramenta em teoria de controle é o uso de controladores não lineares baseados em energia (energy based control ), desenvolvidos através da teoria de Lyapunov (Fauske, 2002 e Khalil, 2000). Fauske (2002) sugere diferentes controladores baseados em energia para satélites equipados apenas com bobinas magnéticas; as leis de controle avaliadas para o modo de estabilização são: 1) realimentação de velocidade angular (angular velocity feedback ); 2) realimentação de atitude (attitude feedback ). Overby (2004) realiza um estudo comparativo dos controladores 1 e 2 com o controlador obtido pela teoria linear ou seja usando a metodologia LQR. Essas estratégias de controle são aplicadas ao satélite Norueguês NCUBE. Hegrenes (2004) realiza um estudo comparativo entre um controlador usando o método LQR e um controlador PD, para o satélite NCUBE, no seu modo de estabilização. A teoria de Lyapunov é desenvolvida em detalhes nas referências Khalil (2000) e 46 Overby (2004). Para o modo de operação de descapotamento/detumble, a lei de controle é sugerida por Wisniewski (1996) e é baseada na teoria de Lypunov. Essa lei de controle também é conhecida como Bdot (Silani e Lovera, 2003). Wells e Jeans (2002) utilizam essa estratégia de controle para o satélite Canadense CanX1. Nesse trabalho a estratégia de controle de Wisniewski (1996) é usada para o modo de operação de detumble. Portanto nos dois sistemas de controle de atitude (ACS(s)), o detumble é feito com o uso de bobinas magnéticas. 3.3 Literatura do INPE No acervo do INPE pesquisado, destacam-se o trabalho de Souza (1981), que apresenta o estudo e desenvolvimento de um sistema de controle ativo de atitude em três eixos para satélites artificiais, usando atuadores pneumáticos a gás frio e rodas de reação. O estudo é feito através de simulação digital; Souza (1987) apresenta o estudo de um sistema de controle de atitude em três eixos utilizando rodas de reação. O autor ainda apresenta o modelo dinâmico do satélite equipado com rodas, utilizando as equações de corpo rı́gido ou equações de Euler, assim como desenvolvido nesse trabalho; Carrara (1982) faz o estudo das principais forças e torques atuantes em satélites artificiais, como gradiente de gravidade, pressão de radiação (solar e terrestre) e arrasto atmosférico, desenvolvendo um programa/software para o computo numérico dessas influências na atitude e órbita do satélite, a partir da sua geometria. Souza (1987) apresenta o estudo e o projeto experimental de um atuador do tipo roda de reação. Sua aplicação está no controle fino da atitude de satélites artificiais atráves da troca de quantidade de movimento angular entre a roda e o satélite, realiza, ainda testes de desempenho e confrontação com os resultados obtidos em simulação digital. Trivelato (1988) apresenta o desenvolvimento de controladores digitais de torques de rodas de reação, usadas no controle de atitude de satélites artificiais, o autor realiza simulações digitais e simulação com hardware na malha de controle, qualificando o controlador da roda para seu uso em sistemas de controle de atitude. O estudo detalhado de um sistema de controle por referência para operação de rodas é apresentado em Trivelato (1988). O projeto de controle de atitude de satélites estabilizados em três eixos utilizando a 47 metodologia LQR e LQG é apresentada em Moscati (1992) mostrando ser uma opção viável para o sistema de controle de atitude, atravês dos resultados obtidos. Moscati (1992) utiliza o modelo linear em suas simulações. Outros trabalhos se destacam no estudo de atitude de satélite que incluem apêndices flexı́veis usando metodologias modernas e distintas de controle (LQG/LTR (Loop Transfer Recovery) e H-infinito), Roma (1991) desenvolve um software de um manipulador simbólico para obtenção das equações do movimento de atitude. Flora () realiza o projeto de um sistema de controle de atitude em três eixos usando os métodos LQG/LTR e H-infinito. Souza (1987) apresenta uma extensão na teoria do regulador linear quadrático (LQR) para o caso de uma dinâmica não linear, sugerindo uma nova lei de controle de atitude para satélites artificiais. A partir da revisão da literatura do INPE nota-se uma grande diversidade e qualidade nos trabalhos voltados para o estudo de atitude em três eixos, e implementações de metodologias de controle. Entretanto, os trabalhos revistos tratam de técnicas que empregam apenas rodas como atuadores, incluindo até mesmo o projeto da mesma. As análises e simulações se restringem a um único modo de operação, o de estabilização, envolvendo apenas manobras de pequenos ângulos. O presente trabalho que inclue o estudo de atitude empregando rodas, especificadas para o satélite EQUARS, é usado para outros modos de operação, estabilização e aquisição. O modelo do satélite equipado com atuadores magnéticos também é desenvolvido. Esse trabalho aborda o problema do uso de rodas usando as metodologias usadas em Souza (1987), Flora () e Moscati (1992) para a obtenção da lei de controle. Na simulação digital foi adotada o modelo completo da dinâmica (não linear e acoplada), assim como feito por Souza (1987). Na literatura disponı́vel do INPE não é encotrada nenhuma referência que trata da utilização de bobinas magnéticas como atuadores para a concepção do sistema de controle de atitude para satélites estabilizados em três eixos. Constitue, portanto, uma importante contribuição e desafio desse trabalho, o estudo da estabilização em três eixos de satélites artificiais através de atuadores magnéticos. Cabe salientar que esse estudo foi motivado pelos projetos dos satélites Brasileiros EQUARS e a Plataforma Multi-Missão (PMM) que estão sendo desenvolvidos pelo INPE. Esses satélites empregam como atuadores rodas e bobinas magnéticas para o controle de atitude. Como contribuição desse estudo desenvolveu-se um toolbox em SIMULINK, voltados para dinâmica e controle de atitude. 48 Neste trabalho as simulações digitais são feitas no ambiente MATLAB e SIMULINK. Nesses recursos computacionais são encontrados nas referências Hanselman e Benjamin (1994), Levine e Leonard (1995), Hanselman e Littlefield (2003) e Bishop (1997) e na documentação disponibilizada pela Mathworks. Levine e Leonard (1995). Ogata (1998) e Bishop (1997) apresentam uma introdução do uso do MATLAB em engenharia de controle e automação. Em Hanselman e Benjamin (1994) e Wie (1998) é apresentado a implementação das metodologias LQR e LQG através das functions do control toolbox do MATLAB. Hanselman e Littlefield (2003) mostra como usar os recursos do MATLAB e apresentam todas as caracterı́sticas e benefı́cios da programação em MATLAB. Hanselman e Littlefield (2003) mostram como realizar o desevolvimento de um toolbox, usando SIMULINK. Para a geração do software de atitude em SIMULINK, foram usados além desses recursos a documentação disponibilizada pela Mathworks (http://www.mathworks.com/support/). 49 50 CAPÍTULO 4 DEFINIÇÕES DE NOTAÇÕES A posição e a atitude de um satélite podem ser representadas de várias formas (Wertz, 1978, Larson e Wertz, 1992). Essa Seção introduz definições e notações necessárias para o desenvolvimento do trabalho proposto. 4.1 Sistemas de Referência Os sistemas de referência utilizados nesse trabalho são definidos a seguir, a Figura (4.1) ilustra os referenciais definidos. 4.1.1 Referencial Inercial O referencial ECI, no qual as leis de Newton do movimento são formuladas, é considerado um referencial inercial (Fauske, 2002). O referencial inercial (X, Y, Z) definido tem origem no centro da Terra. O eixo Z aponta na direção do polo norte geográfico, o eixo X aponta na direção do ponto vernal γ. O eixo Y é encontrado usando a regra da mão direita, completando o sistema dextrogiro. Nota: De fato, o referencial acima é aproximadamente inercial, devido ao movimento da Terra; para as escalas de tempo adotadas, entretanto, o referencial ECI constitui uma boa aproximação de um referencial inercial (Chobotov, 1996). 4.1.2 Referencial Orbital O referencial orbital (OF) definido por (xo , yo , zo ) é um sistema de coordenadas com origem no centro de massa do satélite. O eixo zo aponta na direção do centro da Terra (direção Nadir) e o eixo yo aponta na direção normal a órbita. O eixo xo é obtido pela regra da mão direita, e coincide com a direção do vetor velocidade orbital linear, para uma órbita circular (Wertz, 1978). Nota: Esse referencial também é referido como Vertical Local e Horizontal Local (VLHL). 51 4.1.3 Referencial do Satélite O referencial do corpo (BF), ou do satélite definido por (x, y, z) é um sistema de coordenadas com origem no centro de massa do satélite. Os eixos são escolhidos como sendo coincidentes com os eixos dos momentos principais de inércia. Para estudos de satélite estabilizados em três eixos, Terra-apontado, é prático definir os eixos de roll, pitch e yaw como sendo (Wie, 1998, Moscati, 1992) • eixo de roll x, nominalmente alinhado com xo • eixo de pitch y, nominalmente alinhado com yo • eixo de yaw z, nominalmente alinhado com zo Figura 4.1- Sistemas de referência, inercial (ECI), orbital (OF) e do satélite (BF). 52 4.2 Representação da Atitude A atitude em três eixos é convenientemente representada através de uma transformação de coordenadas, na qual transforma-se um conjunto de coordenadas no espaço inercial em um conjunto de coordenadas fixadas ao satélite (Wertz, 1978). Existem parametrizações/representações alternativas para estas transformações. A principais representações são: • Matriz de rotação ou matriz de atitude; • Ângulos de Euler; • Quaternions ou parâmetros de Euler; • Ângulo eixo equivalente; Cada um desses conjuntos de parametrizações apresentam vantagens e desvantagens. Em Wertz (1978), Hughes (1986), Fauske (2003) e Wie (1998) são discutidas pormenorizadamente cada uma dessas parametrizações. 4.2.1 Matriz de Rotação A matriz de rotação, também chamada de matriz de cossenos diretores, tem as interpretações (Fauske, 2003): • Descreve a orientação mútua entre dois referênciais, onde cada vetor coluna são os cossenos dos ângulos entre os dois referenciais; • Transforma um vetor representado em um referencial para outro; A matriz de rotação R de um referencial a para um referencial b denotada por Rba é um elemento do conjunto SO (3), definido como SO (3) = R | R ∈ <3×3 , RT R = I e detR = 1 , 53 (4.1) onde I é a matriz identidade 3 × 3. A notação seguinte é usada para transformar um vetor r de um referencial para outro: f rom rto = Rto f rom r (4.2) onde o ı́ndice superior descreve em qual referencial o vetor está expresso. Devido a propriedade de ortogonalidade, RT R = I, pode-se mostrar que a derivada no tempo da matriz de rotação é dada por Ṙab = S(ω aab )Rba (4.3) onde a notação usada ω aab representa a velocidade angular do referencial b em relação ao referencial a, expressa no referencial a. Essa notação é usada em todo o trabalho. A velocidade angular tem a propriedade ω aab = −ω aba (Fauske, 2003), e S(ω) é o operador anti-simétrico (skew-symmetric): 0 −ωz ωy S (ω) = ωz 0 −ωx , −ωy ωx 0 ωx ω = ωy ωz (4.4) reescrevemos (4.3), temos Ṙab = S(ω aab )Rba = −S(ω aba )Rba Denotaremos a matriz de rotação por 54 (4.5) Rba = h cb2 cb1 cb3 i (4.6) T onde cbi = cbix cbiy cbiz (i = 1, 2, 3) são matrizes coluna, tendo como componentes os cossenos diretores dos eixos de a em relação aos eixos do referencial b. Daqui para frente o ı́ndice b será usado para representar o referencial do corpo e o ı́ndice o será usado para representar o referencial orbital. A ausência de ı́ndice representa o referencial inercial. 4.2.2 Parâmetros de Euler Os parâmetros de Euler, são também chamados quaternions (Fauske, 2003), e são uma representação atrativa devido a parametrização sem singularidades e equações diferenciais da cinemática serem lineares (Fauske, 2003, Wie, 1998). A representação em quaternion requer menos tempo computacional que a representação por ângulos de Euler (Arantes e Fonseca, 2004a), e portanto é usado em aplicações onde os recursos computacionais são limitados (Wertz, 1978). Um quaternion q é definido como uma grandeza hiper-imaginária composta por uma parte real η e um vetor dada por η = cos Φ , 2 = λ sin Φ 2 (4.7) que, representa uma rotação de Φ em torno do vetor unitário λ. Um quaternion satisfaz o vı́nculo qT q = 1, ou seja η 2 + 21 + 22 + 23 = 1 (4.8) As equações diferenciais da cinemática são dadas por (Wertz, 1978, Wie e Arapostathis, 1989, Wie, 1985): 55 1 q̇ = Ωq 2 (4.9) 0 ωz −ωy ωx −ωz 0 ω ω x y Ω= ω −ω 0 ω x z y −ωx −ωy −ωz 0 (4.10) onde A matriz de atitude (4.1) obtida a partir dos quaternions é dada por (Wie, 1998) Rba = η 2 − qT q I + 2qqT − 2ηQ () (4.11) onde Q () é o operador anti-simétrico dado por 0 −3 2 Q () = 3 0 −1 −2 1 0 (4.12) Dada a matriz de atitude (4.1), podemos determinar η e pelas relações 1 η = (1 + c11 + c22 + c33 ) 2 para 0 ≤ c23 − c32 1 = c31 − c13 4η c12 − c21 Se η = 0 escolhe-se outra sequência de rotações. 56 se η 6= 0 Φ ≤π 2 (4.13) (4.14) 4.2.3 Ângulos de Euler A orientação de um corpo também pode ser descrita por três ângulos (três parâmetros independentes) denominados ângulos de Euler. A Figura (4.2) ilustra a seqüência de rotações necessárias para levar um referencial X̄ ≡ (X, Y, Z) a outro x̄ ≡ (x, y, z), que são listadas como. 0 0 0 • Rotação em torno do eixo Z de um ângulo φ que leva ao sistema ξ , θ , ζ . 0 • Rotação em torno do eixo ξ de um ângulo θ que leva ao sistema ξ, θ, ζ. • Rotação em torno do eixo ζ de um ângulo ψ que leva ao sistema x, y, z. Para essa seqüência de rotações 3 − 1 − 3 a matriz de atitude é dada por (Kaplan, 1976) cos ψ cos φ − sin ψ cos θ sin φ cos ψ sin φ + sin ψ cos θ cos φ sin ψ sin θ RX x = − cos ψ cos φ − cos ψ cos θ sin φ − sin ψ sin φ + cos ψ cos θ cos φ cos ψ sin θ sin θ sin φ − sin θ cos φ cos θ Figura 4.2- Construção dos ângulos de Euler. Em problemas de estabilização de atitude em três eixos é comum definir: • ângulo de roll (φ) é o ângulo de rotação em torno do eixo de roll 57 • ângulo de pitch (θ) é o ângulo de rotação em torno do eixo de pitch • ângulo de yaw (ψ) é o ângulo de rotação em torno do eixo de yaw Note que os ângulos são definidos para rotações em torno de eixos distintos. Usados comumente para relacionar o referêncial orbital (OF) ou LVHL com o referêncial do satélite (BF), adequado para tais aplicações (Moscati, 1992). Existem doze (12) conjuntos possı́veis de ângulos de Euler para descrever um referencial em relação a outro (Wertz, 1978). Esses conjuntos/seqüências são dividas em dois tipos: Tipo 1 (anti-simétrica): Nesse caso as rotações são feitas, sucessivamente, em cada um dos três eixos. Esse tipo apresenta uma singularidade em θ = ± π2 . As seqüências são 1 − 2 − 3, 2 − 1 − 3, 3 − 2 − 1, 2 − 3 − 1, 3 − 1 − 2, 1 − 3 − 2. Tipo 2 (simétrica): Nesse caso a primeira e terceira rotação são feitas sobre o mesmo eixo. Esse tipo apresenta uma singularidade em θ = π e θ = 0. As seqüências são 3 − 1 − 3, 2 − 1 − 2, 1 − 2 − 1, 2 − 3 − 2, 1 − 3 − 1, 3 − 2 − 3. As matrizes de rotação e as equações cinemáticas para cada uma dessas seqüências de rotação (simétricas e anti-simétricas) são encontradas nas referência Wertz (1978) e Hughes (1986). Apesar das equações da cinemática apresentarem singularidades do tipo 1 e tipo 2, os ângulos de Euler tem um clara interpretação fı́sica e são utilizadas na entrada e saı́da das simulações. Dada a matriz de rotação (4.11) escrita em termos de quaternions, os ângulos de roll (φ), pitch (θ) e yaw (ψ) para a seqüência de rotação 3 − 2 − 1, ou seja, primeira rotação em yaw, segunda rotação em pitch e terceira rotação em roll, são dados por (Wie, 1998, Hughes, 1986): φ = arctan c23 c33 0 ≤ φ ≤ 2π π π θ = arcsin (−c13 ) − ≤ θ ≤ 2 2 c12 ψ = arctan 0 ≤ ψ ≤ 2π c11 58 (4.15) (4.16) (4.17) Nota: A seqüencia de rotações escolhida para ser utilizada nesse trabalho é 3−2−1. 4.2.4 Ângulo Eixo Equivalente O ângulo eixo equivalente é definido a partir do teorema de Euler que diz: um corpo rı́gido pode atingir qualquer orientação/atitude a partir de uma orientação arbitrária (de referência), através de uma rotação Φ sobre um eixo unitário fixo λ em relação ao corpo e a referência. Esse eixo é chamado de eixo equivalente ou eixo de Euler. Definindo λ1 λ = λ2 λ3 e 0 −λ3 λ2 Λ = λ3 0 −λ1 −λ2 λ2 0 (4.18) Pode-se mostrar que a matriz de atitude é dada por (Wie, 1998, Fauske, 2002): Rba = cos ΦI + (1 − cos Φ) λλT − sin ΦΛ (4.19) onde I é a matriz identidade. A expressão (4.19) é também conhecida como fórmula de Rodrigues. Dada a matriz de atitude (4.1) Φ pode ser obtido como cos Φ = 1 (c11 + c22 + c33 ) 2 e pode se mostrar (Wie, 1998) que λ é dado por 59 (4.20) c23 − c32 λ1 1 λ = λ2 = c31 − c13 2 sin Φ c12 − c21 λ3 para Φ 6= 0, ±π, ±2π, . . . (4.21) Usando a parametrização (4.19) e a identidade trigonométrica (4.7), a matriz de rotação pode ser escrita em função dos quaternions como Rba = 1 + 2ηS () + 2S2 () 4.2.5 (4.22) Rotações Infinitesimais No modo de estabilização durante a manutenção da atitude nominal, em geral, são realizadas manobras de pequenos ângulos para correções da atitude, justificando aproximar a matriz de atitude por (Moscati, 1992) 1 ψ −θ Rbo = −ψ 1 φ θ −φ 1 (4.23) que é a matriz de transformação do sistema orbital (VLHL) para o sistema do sátelite. Os ângulos φ, θ, ψ são tomados no sentido anti-horário. Para rotações infinitesimais a seqüência de rotações não é importante (Wie, 1998, Moscati, 1992). A representação em quaternions, para pequenos ângulos pode ser aproximado por η∼ =1 φ 1 ∼ = 2 θ 2 ∼ = 2 ψ 3 ∼ = 2 60 (4.24) (4.25) (4.26) (4.27) sendo necessário a normalização dos quaternions para satisfazer a relação η 2 + 21 + 22 + 23 = 1 4.3 (4.28) Sumário de Notações As principais notações usadas nesse trabalho são: Velocidade angular do referencial b em relação ao referencial e exω neb presso/escrito no referencial n ra Vetor r expresso no referencial a a Matriz de rotação do referencial b para o referencial a. ra = Rab rb Rb S (k) Operador anti-simétrico r̂ Versor ou vetor unitário 61 62 CAPÍTULO 5 FORMULAÇÃO DO PROBLEMA Neste Capı́tulo o problema objeto desta dissertação é formulado. A formulação aqui apresentada envolve: • os modelos matemáticos da cinemática e da dinâmica de um satélite contendo 3 rodas de reação e 3 bobinas magnéticas; • a modelagem do gradiente de gravidade e dos torques magnéticos; • a modelagem do campo geomagnético, modelo IGRF; • o tratamento das equações diferenciais, incluindo a linearização dos modelos matemáticos, para serem utilizadas no projeto de controle linear nos modos de controle de estabilização e aquisição; • uma discussão do emprego de atuadores usados neste trabalho com enfoque nas suas aplicações, factibilidade, vantagens e desvantagens; • uma discussão dos torques ambientais incluindo o modelo de pior caso. 5.1 Modelagem Matemática: Cinemática e Dinâmica Esta sessão apresenta e discute os modelos matemáticos da cinemática e da dinâmica de um satélite equipado com três rodas de reação e bobinas magnéticas. As equações da cinemática são obtidas a partir da representação da atitude em ângulos de Euler e quaternions. Essas parametrizações são amplamente utilizadas na área de dinâmica de atitude uma vez que permitem representar a dinâmica e a atitude dos veı́culos espaciais em diferentes sistemas de coordenadas. As equações da cinemática não envolvem as forças ou torques associados ao movimento. Elas descrevem a velocidade do veı́culo espacial em termos de sua orientação em relação a um ou mais sistemas de coordenadas. A orientação do veı́culo em relação a um referencial conhecido é chamada atitude. A atitude pode ser representada por diferentes conjuntos de parâmetros, como é mostrado no Capı́tulo (4). 63 5.1.1 Equações da Cinemática A cinemática descreve a orientação de um veı́culo em relação a um sistema de eixos conhecido. O modelo matemático da cinemática é escrito na forma de um sistema de equações diferenciais de primeira ordem que, em conjunção com as equações de Euler, permite escrever as equações do movimento na forma de estado e, por integração numérica, obter os ângulos de atitude e sua derivada temporal em função do tempo. Em geral as equações da cinemática podem ser representadas por quaternions, matriz de rotação, ângulos de Euler, etc. Cada uma das formas de representar a atitude tem vantagens e desvantagens. A representação por quaternions tem a vantagem de não conter funções trigonométricas e nem conter singularidades, que são tı́picas quando se representa a cinemática em ângulos de Euler. Wertz (1978), Kaplan (1976) e Hughes (1986) apresentam as equações da cinemática em quaternions na forma 1 T b ω ob 2 1 b 1 ˙ = ηω ob − ω bob × 2 2 η̇ = (5.1) (5.2) Usando o operador anti-simétrico S (ω), a Equação (5.2) pode ser reescrita como ˙ = 1 [ηI + S ()] ω bob 2 (5.3) Os parâmetros de Euler ou quaternions dados nas equações (5.2) são apresentados na Seção (4.2). ω bob é o vetor velocidade angular do satélite escrito no referencial do corpo, referido ao sistema orbital. As equações da cinemática podem também ser escritas em função dos ângulos de Euler. A forma do sistema de equações da cinemática em termos dos ângulos de Euler depende da seqüência de rotações escolhidas para se passar de um sistema de referencia para outro. A forma mostrada a seguir pode ser vista em Wertz (1978) e Wie e Arapostathis (1989) e é obtida para a seqüência de rotações 3 − 2 − 1. 64 φ̇ cos θ sin φ sin θ cos φ sin θ sin ψ 1 ωo cos φ cos θ − sin φ cos θ ω bib + θ̇ = 0 cos θ cos ψ cos θ cos θ ψ̇ 0 sin φ cos φ sin θ sin ψ (5.4) onde ωo é a velocidade orbital, que para uma órbita circular é constante. A Figura (5.1) mostra a seqüencia de rotações associadas à Equação (5.4) para representar o sistema do corpo em relação ao sistema orbital, através de rotações subseqüentes dos ângulos ψ, θ e φ. Figura 5.1- Seqüência de rotações 3(ψ) − 2(θ) − 1(φ) . Os eixos x, y e z e os demais sistemas de eixos formados pelas rotações são ortogonais. Entretanto, as rotações em torno dos eixos que são escritas as velocidades angulares não são ortogonais . Por esta razão tem-se seno ou co-seno no denominador (dependendo da seqüência de rotações escolhida) quando se representa o vetor [φ̇ θ̇ ψ̇]T em função dos respectivos ângulos e das componentes do vetor velocidade angular ω bob . Tal caracterı́stica explica a singularidade sempre presente na representação da atitude via ângulos de Euler, uma vez que para cos ±n π2 , n = 1, 3, 5, · · · e sin (±nπ), n = 0, 1, 2, · · · tem-se denominador nulo. A Equação (5.4) ilustra o caso do co-seno. 65 5.1.2 Equações da Dinâmica O modelo real do satélite é aproximado por um modelo fı́sico constituı́do de um corpo rı́gido principal (main bus) contendo internamente 3 rodas de reação, como mostrado na Figura (5.2), também representadas por corpos rı́gidos. O sistema, entretanto, não constitui de fato um corpo rı́gido uma vez que contém partes moveis internas, representadas pelas rodas de reação. O modelo matemático da dinâmica é descrito pelas equações de Euler, expressas no sistema BF (x, y, z) do corpo, referida ao sistema inercial. A origem do sistema de eixos do veı́culo é definida no centro de massa, CM , do satélite. Os eixos x, y e z são coincidentes com os eixos principais de inércia da espaçonave. A Figura ilustra o conceito do satélite objeto da modelagem matemática bem como os sistemas de referências orbital OF(xo , yo , zo ) e do corpo BF (x, y, z). O referencial orbital foi apresentado e comentado na Seção (4.1.2). As equações de Euler utilizadas aqui para representar a dinâmica do satélite, são obtidas a partir da derivada temporal da quantidade de movimento angular, como é mostrado a seguir. Figura 5.2- Ilustração do satélite com rodas de reação, bobinas e os referenciais orbital OF (xo , yo , zo ) e do corpo BF (x, y, z). Seja Lb a quantidade de movimento angular do satélite, expressa no sistema do 66 corpo. O torque externo τext é dado pela taxa de variação no tempo da quantidade de movimento angular. Matematicamente: b dL dLb = + ω bib × Lb = τ ext dt dt b (5.5) onde dLb = dt L̇b = Js ω˙bib + Jw ω˙bib + ω˙bbw (5.6) b onde Js é a matriz de inércia do satélite (sem considerar as rodas), Jw é a matriz de inércia das rodas, ω bbw é a velocidade angular das rodas relativa ao satélite e ω bib é a velocidade angular do satélite em relação ao sistema inercial expresso no sistema do satélite. Lb pode ser escrito com a soma dos vetores quantidades de movimento angular do corpo principal do veı́culo (main bus) mais o vetor quantidade de movimento angular das rodas de reação. Matematicamente Lb = hs + lw (5.7) temos que hs = Js ω bib é a quantidade de movimento angular do satélite e lw = Jw ω bib + ω bbw é a quantidade de movimento angular das rodas. Substituindo na Equação (5.5) Escrevendo Lb dessa forma pode-se reescrever a Equação (5.5) na forma: dLb b = Js ω̇ib + Jw ω̇ bib + ω̇ bbw + ω bib × Js ω bib + Jw ω bib + ω bbw = τ ext dt 67 (5.8) A matriz de inércia Js do corpo principal do satélite (main bus) é dada por Jx −Jxy −Jxz Js = −Jxy Jy −Jyz −Jxz −Jyz Jz Z Jx = 2 y +z B Jxy Z Z x2 + y 2 dm dm, Jy = x + z dm, Jz = B B Z Z Z = (xy) dm, Jxz = (xz) dm, Jyz = (yz) dm 2 2 B 2 B (5.9) (5.10) B Como considerou-se o sistema BF (x, y, z) coincidente com os eixos principais de inércia, essas matrizes são diagonais e os elementos Jx , Jy e Jz são os momentos principais de inércia. Matematicamente: Jx 0 0 Js = 0 J y 0 0 0 Jz (5.11) Jw é matriz de inércia associada às três rodas de reação. Como cada rotor tem seu eixo principal de inércia alinhado com cada um dos eixos principais de inércia do veı́culo, a matriz de inércia Jw , associada aos três rotores, é também diagonal e os elementos da diagonal são os momentos de inércia axiais associados ao eixo de rotação de cada roda de reação. Matematicamente: Jwx 0 0 Jw = 0 Jwy 0 0 0 Jwz (5.12) Reescrevendo a Equação (5.8) em termos das velocidade angulares e sua taxa de variação no tempo na forma: 68 dLb b = (Js + Jw ) ω̇ib + Jw ω̇ bbw + ω bib × (Js + Jw ) ω bib + Jw ω bbw ) = τ ext dt (5.13) A matriz de inércia total do veı́culo é J = Js + Jw . Matematicamente Jx + Jwx 0 0 J = Js + Jw = 0 Jy + Jwy 0 0 0 Jz + Jwz (5.14) Substituindo a Equação (5.14) na Equação (5.13) e usando o operador anti-simétrico, reescrevemos a Equação (5.13) b Jω̇ib + Jw ω̇ bbw + S ω bib Jω bib + Jw ω bbw = τ ext (5.15) Pode-se identificar a expressão correspondente ao torque de controle gerado pelas rodas de reação: Jw ω̇ bbw + S(ω bib )Jw ω bbw = τ bw (5.16) ou seja, o torque das rodas de reação aplicados no satélite para fins de controle e estabilização. Usando esta definição pode-se reescrever a expressão do torque (5.15) como: b Jω̇ib + S ω bib Jω bib = τ ext − τ bw (5.17) O torque externo τ ext representa o somatório de todos os torques externos que atuam no veı́culo espacial. Neste trabalho considerou-se apenas o torque devido ao gradiente de gravidade e o torque de controle gerado pela interação das bobinas magnéticas com o campo magnético da Terra. Considerou-se 3 bobinas, dispostas 69 segundo os eixos x, y e z, como mostrado na Figura (5.3). O torque externo pode, então, ser escrito como na forma: Figura 5.3- Configuração das bobinas magnéticas no satélite. τ ext = τ bd + τ bm (5.18) onde τ bd representa todos os torques de perturbação externa (gradiente de gravidade, pressão de radiação, arrasto atmosféricos, etc) e τ bm representa o torque gerado pelas bobinas, dado por (Fauske, 2002) τ bm = mb × Bb (5.19) onde mb é o momento dipolo magnético gerado pelas bobinas, cujas caracterı́sticas depende de projeto (magnitude e direção do torque) e pode ser representado na forma N x i x Ax mx b b b b m = mx + my + mz = Ny iy Ay = my , N z iz Az mz (5.20) Nessa Equação x, y e z são os eixos do sistema BF (x, y, z) e N , i e A são o número de espiras, a corrente elétrica e a área da espira (parâmetros de projeto), respectivamente, para as três bobinas, dispostas segundo os eixos x, y e z. b B = h Bxb Byb Bzb iT é o vetor campo magnético local, escrito no referencial do 70 satélite. Portanto o modelo matemático que representa o satélite equipado com bobinas é dado por: b Jω̇ib + S ω bib Jω bib = τ bd + τ bm 5.1.3 (5.21) Torque devido ao Gradiente de Gravidade O torque de gradiente de gravidade é causado devido a variação da força gravitacional ao longo do corpo do satélite. O torque depende da altitude, da geometria e da atitude do satélite (Carrara, 1982). Esse efeito é usado para estabilização passiva de satélites O gradiente de gravidade alinha o eixo de menor momento de inércia com a vertical local (Fauske, 2002 e Arantes e Fonseca, 2004a). Esse fato requer que um dos momentos principais de inércia do satélite seja muito menor dos que os outros. Em geral essa configuração é obtida através de mastros extensı́veis colocados no satélite. Uma vez em órbita o mastro é estendido segundo o eixo conveniente a ser alinhado com a vertical local. Exemplos de satélites que usam a estabilização por gradiente de gravidade são DODGE, GEOS-3 e RAE-2. Como deduzido no apêndice A, o torque devido ao gradiente de gravidade, expresso no referencial do satélite, é dado por τ bg = 3ωo2 cb3 × Jcb3 (5.22) onde cb3 é o vetor coluna com os cosenos diretores da vertical local ou eixo zo nas direções dos eixos x, y e z e ωo é a velocidade orbital. 5.2 Campo Magnético Terrestre O campo magnético terrestre é calculado como o gradiente de uma função potencial escalar, que é modelada por um série de harmônicos esféricos, dada por (Fauske, 71 2002) V (r, θ, φ) = a n ∞ X a n+1 X n=1 r m [gnm cos mφ + hm n sin mφ] Pn (θ) (5.23) m=0 A função potencial V é uma função da posição do satélite em coordenadas esfêricas n no sistema ECI definido em 4.1.1. Os coeficientes gm e hnm são chamados coeficientes de Gauss, que são constantes e obtidos empiricamente. A cada cinco anos a IAGA (Association of Geomagnetism and Aeronomy) publica um novo conjunto de coeficientes que constituem o modelo IGRF (International Geomagnetic Reference Field ). Discussões e detalhes desse modelo são encontrados em Macmillan e Quinn (2000). O modelo (código fonte nas linguagens C e FORTRAN) é fornecido em Mclean, . Nesse trabalho a implementação do modelo IGRF é feita no ambiente MATLAB e SIMULINK. A Figura (5.4) mostra as componentes do vetor campo magnético local T Bo = Bxo Byo Bzo , expresso no referencial da órbita (OF) ou VLHL, as componentes do campo magnético são dadas para três revoluções do satélite em torno da Terra. Figura 5.4- Campo magnético local Bo usando o modelo IGRF 2000. A propagação da órbita é descrita no apêndice B. Foi usada uma órbita circular kepleriana. A órbita propagada para o cálculo do campo magnético é a órbita nominal 72 do satélite EQUARS. Os elementos orbitais são mostrados na Tabela (5.1). Tabela 5.1- Elementos orbitais do satélite EQUARS. altura (h) 5.3 excentricidade (e) 750 km ∼ =0 ascenção do nodo ascendente (Ω) 30 graus inclinação (i) 20 graus Tratamento das Equações da Dinâmica para os Casos Estudados nesse Trabalho O controle de atitude de satélites artificiais tem inı́cio na separação foguete-satélite. Até então, qualquer que seja o controle, sua atuação é sobre o veı́culo lançador. Por exemplo, o lançador pode ter um sistema de controle que reduza a rotação do último estágio, para nı́veis de projeto, imediatamente antes da separação. Conhecida a atitude final do último estágio do foguete, tem-se uma estimativa razoável da atitude do satélite imediatamente após a separação. Por outro lado, o último estágio do foguete pode garantir apenas a inserção do veı́culo na órbita correta, não fornecendo dados sobre a atitude do satélite imediatamente após a separação. Nesse caso, o sistema de controle do satélite deve ser capaz de controlar o movimento de atitude do mesmo a partir de qualquer configuração de atitude imediatamente após a separação. Portanto, se o lançador não fornecer, a priori, uma atitude final que possa ser estendida ao satélite imediatamente após a separação, o sitema de controle de atitude deverá, de acordo com requisitos de missões: • Reduzir para valores especificados ou eliminar, de acordo também com especificações de projeto, a rotação do satélite imediatamente após a separação. Esta fase constitui um dos modos de estabilização, comumente conhecido como detumble, traduzido aqui como descapotamento. É a operação de eliminar o movimento de rotação arbitrária do satélite imediatamente após a separação. Este aspecto é tratado nesta Seção porque o controle para este modo pode requerer manobras de grandes ângulos e, neste caso, a linearização do sistema de equações diferenciais e as técnicas de controle linear não se aplicam. Portanto neste modo de operação utiliza-se técnicas de controle não linear; 73 • Uma vez que o sistema de controle de atitude faça o detumbling e tenha controle sobre a velocidade angular e atitude do veı́culo, o satélite entra no modo de estabilização. Uma vez reduzida a velocidade angular do satélite imediatamente após a separação, o satélite deve ser manobrado para adquirir a atitude nominal. Nesse caso, a rotação residual pode sobrecarregar o sistema de controle e, portanto, é desejável que, no modo anterior, a velocidade angular residual não seja crı́tica, ou seja, não comprometa a eficiência do sistema de controle. A fase de estabilização pode requerer, também, manobras de grandes ângulos, uma vez que, no modo detumble, o objetivo é a redução da velocidade angular e não manobras angulares para levar o veı́culo para uma configuração nominal em atitude. Novamente, neste caso, deve se ter um certo conhecimento da atitude do satélite na fase final do detumbling para se decidir sobre a aplicação de leis de controle linear ou não linear. • Uma vez adquirida a condição normal de operação o sistema de controle deve manter o veı́culo estável em torno da atitude nominal, de acordo com especificações de projeto para atender as missões. Nesse trabalho, esse processo também faz parte do modo de estabilização. Nesse modo o controle deve fazer sempre manobras de pequenos ângulos (correções) visando manter o veı́culo na atitude nominal especificada. Neste modo de operação as equações do movimento podem ser linearizadas em torno da atitude nominal e as técnicas de controle linear se aplicam. Neste estudo as técnicas de controle linear aplicadas são as do PD, LQR e LQG. Eventualmente, dependendo da missão, pode-se requerer a mudança de uma atitude nominal para outra, arbitrária, por um certo intervalo de tempo. Nesse trabalho esse modo de contigência é chamado modo de aquisição. Quando isto ocorre novamente requer-se manobras de grandes ângulos e posterior reaquisição da nova atitude nominal. A visão do sistema de controle em modos de operação é mais prática sob o ponto de vista de que o projeto de sistemas de controle reais é implementado por modos de operação. São necessárias diferentes estratégias e técnicas de controle para cada modo. 74 5.4 Linearização Para o projeto do controlador linear usando as metodologias LQR e LQG, que são baseadas em uma planta linear, é feita a linearização das equações do movimento em torno de um ponto de operação. As equações linearizadas também podem ser usadas para estudos de estabilização/manutenção da atitude nominal, para os casos onde são realizadas manobras de pequenos ângulos (< 15o ). As equações usadas são escritas em termos dos ângulos de Euler, já que a sigularidade do tipo 1 (±π/2), não acontece para a faixa de manobras realizadas (pequenos ângulos). Outra justificativa é a clara interpretação fı́sica que os ângulos de Euler fornecem. São feitas as linearizações para as versões: • Satélite equipado com rodas, Equacão (5.17); • Satélite equipado com bobinas, Equação (5.21). A Equação da cinemática (5.4) pode ser escrita como (Bryson, 1994) cos θ sin ψ 1 0 − sin θ φ̇ ω bib = 0 cos φ sin φ cos θ θ̇ − ωo sin φ sin θ sin ψ + cos φ cos ψ cos φ sin θ sin ψ − sin φ cos ψ ψ̇ 0 − sin φ cos φ cos θ (5.24) Para pequenos ângulos a expressão (5.24) pode ser aproximada por (Bryson, 1994) ω1 ≡ φ̇ − ωo ψ − ψ̇θ ω2 ≡ θ̇ − ωo − ψ̇φ ω3 ≡ ψ̇ + ωo φ − θ̇φ (5.25) O valor médio, por órbita, dos termos não lineares ψ̇θ, ψ̇φ, θ̇φ são pequenos comparados com as magnitude média dos termos lineares em (5.25). Se φ e θ oscilam em torno de zero ou se a magnitude de ψ̇ e θ̇ são pequenos comparados com a velocidade orbital ωo . Nesses casos a Equação (5.25) pode ser aproximada por 75 φ̇ ≡ ω1 + ωo ψ θ̇ ≡ ω2 + ωo ψ̇ ≡ ω3 − ωo φ Usando a notação ω bib = [ω1 ω2 (5.26) ω3 ]T temos φ̇ −ψ b ω ib = I θ̇ + ωo −1 ψ̇ φ (5.27) onde I é a matriz identidade de ordem 3. 5.4.1 Linearização do Modelo do Satélite Equipado com Rodas Definimos o vetor de variáveis de estado como h iT x = φ θ ψ φ̇ θ̇ ψ̇ (5.28) Usando a notação hw = [hw1 hw2 hw3 ]T onde hwi , (i = 1, 2, 3) são as componentes, ao longo das direções (x, y, z), da quantidade de movimento angular total das rodas. Para a configuração adotada para as rodas, cada componente da quantidade de movimento angular corresponde a uma roda. Substituindo a Equação (5.27) na Equação (5.17) e incluindo o torque de gradiente de gravidade dado pela Equação (5.22), obtemos Jx φ̈ = [4ωo2 (Jz − Jy ) + ωo hw2 ] φ + ωo h3 + [ωo (Jx − Jy + Jz ) + hw2 ] ψ̇ + [ωo (Jy − Jz ) φ + hw3 + (Jz − Jy )] θ̇ − ḣw1 76 (5.29) Jy θ̈ = 3ωo2 (Jz − Jx ) θ h− ωo φh1 − [ωo2 (Jz − Jx ) φ + ωo hw3 ]iψ − [ωo (Jz − Jx ) ψ + hw1 ] ψ̇ − ωo (Jx − Jz ) φ − hw3 + (Jx − Jz ) ψ̇ φ̇ − ḣw2 (5.30) Jz ψ̈ = [ωo2 (Jxh− Jy ) + ωo hw2 ] ψ − ωo hw1 − [ωo (Jxi− Jy + Iz ) + hw2 ] φ̇ + ωo (Jy − Jx ) ψ + hw1 + (Jx − Jy ) φ̇ θ̇ − ḣw3 (5.31) onde Jx , Jy e Jz são os momentos principais de inércia. Transformando as equações (5.29), (5.30) e (5.31), de segunda ordem em equações de primeira ordem e fazendo a expansão em série de Taylor em (x = 0), ou seja, a vertical e horizontal local (VLHL), obtemos x∗ = h∗w = 0 h 0 Ho iT 0 (5.32) Onde o termo Ho é a quantidade de movimento angular na direção do eixo de pitch gerado pelo volante de inércia. Esse procedimento é utilizado para produzir rı́gidez giroscópica durante a fase de manutenção da atitude. Escrevendo na forma de estados as equações linearizadas, temos ẋ = Ax + Bu com 77 (5.33) 0 0 0 A = 4ωo2 (Jz −Jy )−ωo Ho Jx 0 0 0 0 0 0 3ωo2 (Jz −Jx ) Jy 0 0 0 0 0 0 1 0 0 0 0 ωo2 (Jx −Jy )−ωo Ho Jz ωo (−Jx +Jy −Jz )+Ho Jz 0 1 0 0 0 0 0 1 ωo (Jx −Jy +Jz )−Ho Jx 0 0 0 (5.34) e 0 0 0 0 0 0 B= 1 − J 0 x 0 −1 Jy 0 0 0 0 0 0 0 − J1z (5.35) O vetor de controle é dado por iT h u = u1 u 2 u3 (5.36) u1 = ḣw1 − ωo hw3 u2 = ḣw2 u1 = ḣw3 + ωo hw1 (5.37) onde 5.4.2 Linearização do Modelo do Satélite Equipado com Bobinas Definimos o vetor de variáveis de estado e o vetor de controle como 78 h x = φ θ ψ φ̇ θ̇ ψ̇ iT (5.38) e h u = mx my mz iT (5.39) Usando a notação Bb = Bxb Byb Bzb , onde Bib ,(i = x, y, z) são as componentes do vetor campo magnético local expresso no sistema do satélite. Substituindo a Equação (5.27) na Equação (5.21) e incluindo o torque de gradiente de gravidade dado pela Equação (5.22), obtemos Jx φ̈ = h 4ωo2 i (Jz − Jy ) − ωo (Jz − Jy ) θ̇ φ + ωo (Jx − Jy + Jz ) ψ̇ + (Jy − Jz ) θ̇ψ̇ +my Bzb (t) − mz Byb (t) (5.40) i h Jy θ̈ = 3ωo2 (Jx − Jz ) θ − ωo2 (Jz − Jx ) ψ + ωo (Jx − Jz ) φ̇ φ + ωo (Jx − Jz ) ψ ψ̇ + (Jz − Jx ) φ̇ψ̇ + mz Bxb (t) − mx Bzb (t) (5.41) h i h i Jy ψ̈ = ωo2 (Jx − Jy ) + ωo (Jy − Jx ) θ̇ ψ − ωo (Jx − Jy + Jz ) + (Jy − Jx ) θ̇ φ̇ +mx Byb (t) − my Bxb (t) (5.42) Transformando as equações (5.40), (5.41) e (5.42), de segunda ordem, em equações de primeira ordem e fazendo a expansão em série de Taylor em torno de (x = 0), ou seja, a vertical e horizontal local (VLHL), obtemos ẋ = Ax + B (t) u 79 (5.43) com 0 0 0 A = 4ωo2 (Jz −Jy ) Jx 0 0 0 0 0 0 3ωo2 (Jz −Jx ) Jy 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 ωo2 (Jx −Jy ) Jz ωo (−Jx +Jy −Jz ) Jz 0 0 1 ωo (Jx −Jy +Jz ) Jx 0 0 0 (5.44) e 0 0 0 B= 0 b − Bz (t) Jy Byb (t) Jz 0 0 0 Bzb (t) Jx 0 B b (t) − Jyz 0 0 0 Byb (t) − Jx Bxb (t) Jy 0 (5.45) O sistema linear dado pela Equação (5.43) é variante no tempo (LTV), devido ao campo magnético terrestre não ser constante. O campo é aproximadamente periódico. No sistema linearizado, o vetor campo magnético expresso no referencial do satélite (Bb ) coincide com o vetor campo magnético expresso no referencial orbital (Bo ). Portanto a matriz B pode ser escrita como 0 0 0 B= 0 o − Bz (t) Jy Byo (t) Jz 0 0 0 Bzo (t) Jx 0 B o (t) − Jyz 80 0 0 0 o By (t) − Jx Bxo (t) Jy 0 (5.46) As equações (5.43) e (5.33) representam o modelo matemático linearizado do satélite equipado com rodas (de reação e volante de inércia) e bobinas magnéticas, que constituem uma boa aproximação do modelo para manobras de pequenos ângulos (< 15o ). Os projetos do controladores multivariável LQR e LQG foram feitos sobre o modelo linear dado pelas equações (5.43) e (5.33). Entretanto, a simulação digital é feita sobre o modelo completo, não linear, equações (5.4), (5.21) e (5.17). O sistema de controle no modo normal de operação (que é a fase de manutenção da atitude nominal), constitue um problema de regulador, ou seja, a referência é constante e igual a zero (x = 0). O que corresponde ao satélite estar alinhado com a VLHL. O modelo linear do satélite equipado com bobinas, a atitude é descrita em quaternions e é dado em Arantes e Fonseca (2004b) e Overby (2004). No modelo do satélite parametrizado em quaternions o estado é dado por h x = 1 2 3 ˙1 ˙2 ˙3 iT (5.47) Escrevendo o modelo na forma da Equação (5.43), as matrizes A e B são dadas por 0 0 0 A = 4ωo2 (Jz −Jy ) Jx 0 0 0 0 0 0 3ωo2 (Jz −Jx ) Jy 0 0 0 0 0 0 1 0 0 0 0 ωo2 (Jx −Jy ) Jz ωo (−Jx +Jy −Jz ) Jz e 81 0 1 0 0 0 0 0 1 ωo (Jx −Jy +Jz ) Jx 0 0 0 (5.48) 0 0 0 B= 0 o − Bz (t) 2Jy Byo (t) 2Jz 0 0 0 Bzo (t) 2Jx 0 Byo (t) − 2J z 0 0 0 Byo (t) − 2Jx Bxo (t) 2Jy 0 (5.49) Nota: As linearizações das equações da dinâmica foram feitas com o auxı́lio do manipulador simbólico do programa MATHEMATICA e as equações obtidas foram aferidas com as encontradas na literatura. A próxima Seção discute algumas considerações sobre os modelos matemáticos e sistemas de referencias associados, e os atuadores escolhidos para compor os sistemas de controle nesse trabalho. 5.5 Considerações Sobre os Atuadores, os Modelos Matemáticos e Sistemas de Referênica Adotados As considerações sobre os atuadores são mostradas a seguir. 5.5.1 Modelo dos Atuadores Existem vários tipos de atuadores que podem ser usados para o controle de atitude de satélites artificiais. Os atuadores podem ser dividos em três categorias: 1) propulsores; e 2) dispositivos de troca de quantidade de movimento angular, como rodas e 3) atuadores magnéticos. É comum o emprego de mais de um tipo de atuador em satélites, dependendo dos requisitos da missão, e das caracterı́sticas do projeto de controle, Wertz (1978) e Larson e Wertz (1992) apresentam uma descrição do emprego de atuadores, empregados em diferentes missões. Roda de Reação e Volante de Inércia As rodas de reação e volantes de inércia podem ser modelados como motores DC, essa modelagem pode ser encontrada em Bryson (1994), Wilson (2000) e Trivelato 82 (1988). Uma aproximação para a função de transferência desses atuadores é dada em Souza (1981). Define-se (Larson e Wertz, 1992) • Rodas de reação: rodas com quantidade de movimento angular nominal nulo • Volante de inércia: são rodas de reação com quantidade de movimento angular nominal diferente de zero Satélites equipados com volantes de inércia são referidos como satélites com quantidade de movimento angular embarcado (momentum bias system), e satélites equipados com rodas de reação são referidos como satélites com quantidade de movimento angular nulo (zero momentum system). Basicamente, o funcionamento consiste na geração de torques devido a aceleração de uma roda, ligada ao rotor de um motor elétrico, em relação a seu estator que é fixado à estrutura do satélite (Souza, 1987). A dinâmica da roda é função do seu atrito, sua velocidade, e voltagem aplicada pelo motor. Considera-se nesse trabalho que exista um sistema de controle fino para a operação das rodas, que seja eficiente, fornecendo o torque requisitado pelo sistema de controle de atitude. Um estudo detalhado de um sistema de controle por referência para operação de rodas pode ser encontrado em Trivelato (1988). Em Souza (1987) é apresentado o estudo e o projeto de um modelo experimental de uma roda de reação. O uso de volantes de inércia, que apresentam quantidade de movimento angular nominal diferente de zero, fornece rı́gidez giroscópica ao satélite. Para as rodas de reação espera-se que operem em torno de velocidade (relativa) nula ao longo de um perı́odo orbital, adequada para absorver torques externos cı́clicos. A presença de torques seculares, pode levar à velocidade angular máxima da roda, exigindo a sua dessaturação. Isso pode ser feito por bobinas magnéticas e/ou propulsores. Nesse trabalho considerou-se somente rodas de reação combinadas com atuadores magnéticos para compor o sistema de controle. Propulsores não foram considerados. 83 Bobinas Magnéticas As bobinas magnéticas podem ser modeladas como um circuito RL (resistência e indutância) usando as leis de Kirchoff. As bobinas tem uma dinâmica muito rápida comparada com a dinâmica do satélite. Incluindo a dinâmica das bobinas na simulação, o sistema torna-se muito lento, pois envolve duas escalas de tempo distintas, portanto, a dinâmica das bobinas é ignorada na simulação. O erro indroduzido devido a isso é negligenciável (Fauske, 2002). Logo, a dinâmica total do sistema é relativamente lenta. Bobinas magnéticas são comumente usadas em aplicações espaciais desde que os requisitos de tempo não sejam muito exigentes. As manobras de grandes ângulos via bobinas magnéticas podem levar horas, dependendo das especificações da bobina. As bobinas são também muito usadas para a dessaturação de rodas de reação. A saturação das rodas ocorre quando elas atingem sua quantidade máxima de quantidade de movimento especificada. Quando isto ocorre não se consegue mais acelerar as rodas e, portanto, não se consegue mais usá-las para finalidade de torque objetivando o controle da atitude. As bobinas podem também ser usadas para manobras de pequenos ângulos. Uma aplicação muito comum para manobras de grandes ângulos se refere as operações do controle no modo detumble. O estudo de dessaturação das rodas pode ser encontrado nas referências Bryson (1994) e Bang e Choi (2003). Para manobras de pequenos e grandes ângulos, ou seja, estabilização e aquisição, respectivamente, as bobinas são, em geral, usadas para pequenos satélites (< 500kg) (Wisniewski e Blanke, 1999,Wisniewski e Markley, 1999), como é o caso do satélite brasileiro EQUARS (∼ = 150kg). O controle de pequenos satélites usando somente bobinas magnéticas só é possı́vel para órbitas em que a variação do campo magnético seja suficiente para garantir o controle do satélite (Wisniewski, 1997, Musser e Ebert, 1989). O uso das bobinas magnéticas para operações em órbita faz uso da interação entre o momento dipolo magnético gerado pelas bobinas e do campo magnético terrestre. Por isso, para efeitos de simulações para avaliação do desempenho dos torques magnéticos é muito importante um bom modelo do campo magnético terrestre. Em termos de projeto é sempre possı́vel, dentro das limitações de espaço, massa e energia, alterar o momento das bobinas magnéticas, via parâmetros de projeto tais 84 como número de espiras, corrente disponı́vel e área das bobinas. Entretanto não se tem controle do campo geomagnético, que é parte do ambiente espacial. O modelo do campo geomagnético e sua eficiência para os sistemas de controle com bobinas magnéticas depende dos parâmetros orbitais. Portanto, o estudo do controle via bobinas exige conhecimento do modelo do campo bem como do modelo da órbita do satélite em questão. Esse Capı́tulo é concluı́do com a apresentação dos modelos de pior caso para os principais torques ambientais e um modelo de torques de perturbação internos. 5.6 Torques Ambientais Um satélite está sujeito a pequenos mas a persistentes torques externos devido a perturbações ambientais. Sem resistência a essas perturbações o satélite pode perder a atitude nominal. Os principais torques ambientais são: • torque devido ao gradiente de gravidade; • torque aerodinâmico; • torque devido a pressão de radiação (solar e terrestre); • torque devido a interação do dipolo residual com o campo magnético. Essa Seção discute brevemente esses torques ambientais, e apresenta o modelo de pior caso, conveniente para dimensionalização dos atuadores (torque de controle), para satélites em geral. Uma discussão pormenorizada é encontrada em Carrara (1982), Wertz (1978), Hughes (1986) e Fauske (2002). Nesse trabalho somente o torque devido ao gradiente de gravidade é considerado. Esse torque é considerado, nesse estudo, como o mais importante em termos de perturbação, desde que se considere um bom balanceamento magnético. Os outros torques não são considerados. 5.6.1 Torque Devido ao Gradiente de Gravidade Este torque depende da altitude do veı́culo, das suas propriedades de inércia e de sua atitude em relação à vertical local. Uma expressão que leva em conta estes aspectos e 85 que é apropriada para se calcular em primeira aproximação a magnitude do torque do gradiente de gravidade para o pior caso é dada em Larson e Wertz (1992) e reproduzida aqui τgb = 3ωo2 |Jx − Jz | sin 2θ 2 (5.50) onde θ é o ângulo que descreve o deslocamento angular na direção da vertical local ou nadir. O gradiente de gravidade tem a propriedade de alinhar o eixo de menor momento de inércia com a vertical local, em uma configuração de estabilidade chamada estabilização por gradiente de gravidade. Por esta razão satélites estabilizados por gradiente de gravidade requerem o projeto estrutural de tal forma a permitir que o eixo de apontamento para a Terra seja o eixo de menor momento de inércia. Um exemplo interessante refere-se à primeira versão do satélite brasileiro de coleta de dados, SCD-1, que deveria ser estabilizado por gradiente de gravidade. O veı́culo deveria conter um mastro a ser destendido em órbita. O mastro teria um comprimento de 10 m, com uma massa de 3 kg na ponta. Isto aumentaria de aproximadamente dez vezes os momentos de inércia transversais do satélite enquanto manteria o momento de inércia axial (eixo de apontamento para a Terra) aproximadamente igual. Posteriormente a configuração de estabilização foi reconsiderada e o projeto modificado para a estabilização por spin (direção de ω coincidente com eixo axial do veı́culo). Esta configuração de estabilização requer apenas que o eixo de rotação do veı́culo seja o de maior momento de inércia. A mudança diminuiu os riscos associados as operações de abertura e captura do satélite (SCD-1) em torno da vertical local. Normalmente os satélites estabilizados por gradiente de gravidade requerem um mastro a ser aberto em órbita para adequar as propriedades de inércia do satélite a estabilização por gradiente de gravidade. Portanto a complicada dinâmica na fase transitória da abertura dos mastros representam um aspecto delicado para tal tipo de estabilização. Como existem restrição de espaço e volume nos veı́culos lançadores não há como lançar um satélite estabilizado por gradiente de gravidade já com a configuração adequada em termos de propriedades de inércia. A geometria em forma de lápis (pencil-shaped ) que caracteriza os satélites estabilizados por gradiente de gravidade é, de fato, inadequada para os veı́culos lançadores. Daı́ a necessidade de mastros com massa adicional na ponta para dar ao satélite a forma e as propriedades de inércia apropriadas para a estabilização por gradiente de 86 gravidade. Neste trabalho considerou-se apenas o torque do gradiente de gravidade como perturbação ambiental. A justificativa é que assumi-se um satélite de órbita baixa da Terra, LEO (Low Earth Orbit). Assumiu-se também que o satélite objeto de simulações, neste caso o EQUARS, tenha um balanço magnético adequado de tal forma que a perturbação predominante seja associada ao gradiente de gravidade. Entretanto outros torques podem perturbar o movimento do satélite em órbita da Terra. A seguir ao feitas algumas considerações sobre outros torques importantes na análise da dinâmica e controle de atitude de satélites artificiais da Terra. 5.6.2 Torque Aerodinâmico Para satélites de órbitas baixas (< 900km) a densidade do ar é suficiente para influenciar a dinâmica de atitude do satélite. Além da altitude o arrasto atmosférico também depende das dimensões do satélite, da sua geometria e velocidade relativa (Fauske, 2002). O torque aerodinâmico pode ser escrito como 1 τ a = ρV 2 Cd A (uv × scp ) 2 (5.51) onde ρ é a densidade da atmosfera em (kg/m3 ), A em (m2 ) é a área perpendicular a uv , uv é o versor na direção da velocidade, Cd é o coeficiente de arrasto aerodinâmico, V é a velocidade em (m/s) e scp é o vetor distância do centro de massa ao centro de pressão. Uma expressão para o pior caso é dada por Fauske (2002) τa = F (cpa − cg ) F = 0.5 ρCd AV (5.52) 2 onde cg é o centro de gravidade e cpa é o centro de pressão. 87 (5.53) 5.6.3 Torque de Pressão de Radiação Solar A radiação e partı́culas do Sol afetam a dinâmica do satélite. Para baixas órbitas o efeito pode ser negligenciado se comparado com outros torques ambientais (Fauske, 2002). Uma expressão para o pior caso é dada por Larson e Wertz (1992) τs = F (cps − cg ) Fs F = As (1 + q) cos i c (5.54) (5.55) onde Fs é uma constante solar (1367W/m2 ), As é a área superficial irradiada em (m2 ), cg é o centro de gravidade (m), i é o ângulo de incidência do Sol, c é a velocidade da luz na vácuo em (m/s), cps é o centro de pressão solar (m) e q é o fator de reflectância, que varia de 0 a 1. 5.6.4 Torque Devido ao Dipolo Residual Os sistemas eletrônicos no satélite podem gerar um dipolo magnético residual. Esse dipolo residual irá interagir com o campo magnético da Terra. O dipolo residual também mascara as medidas de Bb , feitas pelos magnetômetros. O torque resultante pode ser expresso como (Fauske, 2002) τm = D · B (5.56) onde D é o dipolo magnético residual em (Am2 ) varia de 0.1 à 20Am2 , dependendo do tamanho e desbalanceamento magnético do satélite, e B é o campo magnético da Terra medido em Tesla. Para uma órbita polar B pode ser aproximado por B= 2M Rs2 88 (5.57) onde M = 7.95 · 1015 T · m3 é o momento magnético da Terra e Rs é distância do centro de massa da Terra ao centro de massa do satélite. Nesta Seção faz-se também algumas considerações sobre torques internos ao sistema. Esses torques não alteram a quantidade de movimento angular do sistema, implicando torque externo nulo. Entretanto, os efeitos dos torques internos podem causar problemas no apontamento de antenas, painéis, telescópios (tipo Hubble). Por outro lados efeitos internos dissipadores podem ser úteis, principalmente quando se tratam dos amortecedores de nutação, quase sempre requeridos quando se usa estabilização por spin (SCD-1 e SCD-2). 5.7 Perturbações internas Perturbações internas são torques internos exercidos sobre o corpo do satélite devido a partes móveis. O efeito de torques internos é a dissipação de energia cinética e a redistribuição das componentes da quantidade de movimento angular no sistema do satélite (Wertz, 1978), como acontece na presença de rodas. Entretanto, a presença de perturbações internas não altera a quantidade de movimento angular referido ao espaço inercial. As principais fontes de perturbações internas são: • Movimento de painéis, rodas, mastros, antenas, etc; • slosh devido ao movimento de lı́quidos (combustı́vel); • Deformação na estrutura (snap) devido a dilatação térmica; • Choque de gases no corpo do satélite propelidos pelo sistema de propulsão; • etc. Um exemplo bem conhecido onde a influência de perturbações internas, provocado pelo movimento de antenas, é o do satélite explorer 1, lançado em 1958. O efeito de dissipação de energia, resultou na mudança do eixo de rotação, eixo de menor momento de inércia, para o eixo de maior momento de inércia (Kaplan, 1976). Perturbações internas são difı́ceis de serem modeladas e podem gerar jitter. Uma aproximação para esse efeito é fornecida por Wilson (2000) 89 T = a0 + X ak cos (kωo t + β) + bk sin (kωo t + β) (5.58) k no qual ωo é o velocidade orbital, ao , ak e bk são coeficientes constantes e β a fase. A influência desse efeito é conveniente para estudos de longa duração, particularmente para sistema de controle com quantidade de movimento angular embarcado (Wilson, 2000). 90 CAPÍTULO 6 PROJETO DE CONTROLE Este Capı́tulo é dedicado a apresentar as estratégias de controle usadas para os modos de operação do satélite: 1) detumble, fase em que o satélite é injetado em órbita; 2) estabilização, fase em que o satélite adquire a vertical local e realiza a manutenção da atitude nominal e 3) aquisição, fase em que o satélite faz uma reorientação da atitude. É apresentado o controlador Bdot, os controladores baseados em energia e as metodologias do Regulador Linear Quadrático (LQR) e do Regulador Quadrático Gaussiano (LQG) usados no modo de estabilização. O método LQR é analisado primeiramente, formando a base para um melhor entendimento do método LQG. Essas estratégias são discutidas pormenorizadamente nas referências Dorato e Cerone (1995), Maciejowski (1989), Kwakernaak e Sivan (1972), Moore e Anderson (1990), Brown (1997), Kirk (1970), Fauske (2002) e Overby (2004). 6.1 Modo de Detumble O objetivo do controlador usado no modo de detumble é dissipar a energia cinética do satélite, reduzindo sua velocidade de rotação. 6.1.1 Controlador Bdot O controlador usado nesse trabalho para o detumble é proposto por Wisniewski (1996) e é conhecido como controlador Bdot (Silani e Lovera, 2003). O controlador Bdot ou de Wisniewski (1996) usa a taxa de variação das medidas feitas pelos magnetômetros embarcados no satélite. Sua viabilidade é confirmada pela aplicação no satélite Canadense CanX-1 (Wells e Jeans, 2002). A lei de controle é dada por mb = −k Ḃb − mc onde 91 (6.1) mc = [0 0 mc ] (6.2) no qual k é uma constante escalar maior que zero e mc o momento dipolo nominal da bobina disposta na direção z. A demostração analı́tica de que a lei de controle Bdot dissipa a energia cinética de rotação do satélite, é dada usando a teoria de Lyapunov. Essa demostração é encontrada em Fauske (2002) e Wisniewski (1996). Nos dois procedimentos ou sistemas de controle de atitude propostos, o detumbling é feito com o uso do controlador acima, através dos atuadores magnéticos. 6.2 Modo de Estabilização Para modo de estabilização definido no Capı́tulo 1 as leis de controle são baseadas: • na metodologia LQR e LQG; • nos controladores baseados em energia. Essas estratégias são apresentadas a seguir. 6.2.1 Método LQR O método LQR é baseado na linearização de sistemas dinâmicos, pois a metodologia é formulada para sistemas lineares. O problema de controle ótimo consiste em minimizar uma função custo quadrática e gerar uma matriz de ganhos para realimentação (Dorato e Cerone, 1995 e Maciejowski, 1989). A metodologia LQR é formulada para os dois sistemas dinâmicos envolvidos: 1) satélite equipado com rodas e 2) satélite equipado com bobinas. As equações linearizadas dos modelos, Equação (5.33), invariante no tempo e Equação (5.43) variante no tempo são utilizadas para o projeto de controle. A formulação do problema é feita no domı́nio do tempo e é descrita a seguir. 92 Seja um sistema linear descrito por ẋ = A(t)x + B(t)u (6.3) O problema de otimização consiste em encontrar uma lei de controle linear do tipo u = −Kc (t)x (6.4) que minimize o ı́ndice de desempenho quadrático Z Jp = T T x Qc x + uT Rc u (6.5) 0 onde Qc é uma matriz positiva semi-definida (Qc ≥ 0), Rc positiva definida (Rc > 0), x é o vetor de estado de dimensão n × 1 e u é o vetor de controle de dimensão m × 1. As matrizes Qc e Rc são as ponderações no vetor de estado e no vetor de controle, respectivamente. Nota: para a existência e estabilidade da solução do problema LQR, a condição necessária e suficiente é que o sistema seja completamente controlável (Dorato e Cerone, 1995 e Maciejowski, 1989). As análises de controlabilidade, mostradas em Kwakernaak e Sivan (1972), feitas para os sistema envolvidos. As equações (5.33) e (5.43), garantem a solução do problema LQR. A solução do problema LQR, ou seja, o cálculo do ganho de controle é obtido resolvendo a Equação de Ricatti, dada por (Moore e Anderson, 1990) T Ṗc = −Pc A(t) − A(t)T Pc + Pc B(t)R−1 c B Pc − Qc (6.6) Resolvendo a Equação diferencial matricial de Ricatti obtemos o controlador variante no tempo: 93 u = −R−1 c B(t)Pc (t)x (6.7) Kc = R−1 c B(t)Pc (t) (6.8) Da Equação (6.4), temos A Figura (6.1) mostra a configuração do LQR para o sistema em malha fechada. Figura 6.1- Configuração do controle LQR. O controlador (6.7) exige um grande esforço computacional, devido ao fato de que a cada passo de integração Pc (t) deve ser computado, nos casos em que o sistema for variante no tempo, Equação (5.43). Método LQR Aplicado para o Satélite Equipado com Rodas Devido ao fato do sistema (5.33) ser invariante no tempo e considerando o intervalo de otimização infinito, ou seja, T da Equação (6.5) tendendo ao infinito. A solução do problema LQR é dada para o estado estacionário. Esse problema é referido como problema do LQR estacionário (steady-state LQR) (Dorato e Cerone, 1995). Nessas condições a Equação diferencial matricial de Ricatti torna-se uma Equação matricial algébrica, dada por: T 0 = −Pc A − AT Pc + Pc BR−1 c B Pc − Qc 94 (6.9) e a lei de controle para o sistema dado pela Equação (5.33) resulta em u = −R−1 c BPc x (6.10) onde as matrizes A e B são dadas em (5.34) e (5.35), respectivamente. Método LQR Aplicado para o Satélite Equipado com Bobinas Para a formulação da lei de controle usada no sistema (5.43) utilizando a metodologia LQR, pode ser feita a seguinte aproximação: consideramos o campo magnético apresentado na seção 5.2, periódico, e tomamos um valor médio de Bo . Dessa forma assumimos que o modelo (5.43) seja invariante no tempo. Com essas considerações e tomando (T → ∞) o problema se reduz ao problema LQR estacionário. Temos que a Equação de Ricatti resulta T 0 = −Pc A − AT Pc + Pc BR−1 c B Pc − Qc (6.11) e a lei de controle para o sistema dado pela Equação (5.33) resulta em u = −R−1 c BPc x (6.12) Overby (2004) e Silani e Lovera (2003) apresentam resultados satisfatórios usando o procedimento acima. Entretanto, nesse trabalho, a Equação matricial algébrica de Ricatti é calculada para cada passo de integração, considerando-se a matriz de entrada B(t) variante no tempo, obtida a partir das componentes do campo geomagnético, dada pela Equação (5.49). A solução da Equação algébrica matricial de Ricatti resulta em um ganho variante no tempo para o controlador, dado por u = −R−1 c B(t)Pc (t)x 95 (6.13) ou u = −Kc (t)x (6.14) onde Kc (t) = R−1 c B(t)P(t). As matrices de ponderação Qc e Rc são definidas como: Rc = diag ([r1 , r2 , · · · , rna ]) (6.15) Qc = diag ([q1 , q2 , · · · , qns ]) (6.16) onde na é o número de atuadores no sistema de controle e ns é o número de estados de interesse. Para o EQUARS nos dois procedimentos/modelos temos na = 3 e ns = 6. O desempenho desejado do sistema é obtido pelo ajuste das ponderações. Como sugerido por Overby (2004) e Kristiansen (2000) as ponderações são escolhidas como: qi = 1 (∆xi )2 e ri = 1 (∆ui )2 (6.17) Os valores de ∆ui são baseados no máximo esforço de controle ou valor máximo de operação dos atuadores, ou seja, máximo torque no caso de rodas e máximo momento dipolo no caso de bobinas. Os valores de ∆xi são baseados na faixa/intervalo de operação dos estados. Para os sistemas dinâmicos avaliados as escolhas dos parâmetros usados são baseadas nas especificações do satélite EQUARS, temos ∆xi = 10 graus, ∆ui = 75 mN m (i = 1, 2, 3) e ∆ẋi = 1 graus/s (i = 4, 5, 6) (6.18) 96 para o modelo do satélite equipado com rodas e ∆xi = 10 graus, ∆ui = 1 Am2 (i = 1, 2, 3) e ∆ẋi = 1 graus/s (i = 4, 5, 6) (6.19) para o modelo do satélite equipado com bobinas. 6.2.2 Método LQG Nessa seção o problema do LQG é apresentado. A metodologia será implementada ao sistema (5.43). A solução é obtida no estado estacionário (T → ∞). As medidas são corrompidas por ruı́dos, modelados como distribuições aleatórias gaussianas. É apresentado a estrutura do filtro de Kalman-Bucy e o princı́pio da separação. O problema LQG é apresentado como segue. Dado o sistema da Figura (6.2), onde G é a função de transferência da planta. O problema LQG pode ser colocado como sendo o de calcular uma lei de controle que matenha o sistema estável e minimize um critério de erros quadráticos (Maciejowski, 1989). O sinal de referência r na Figura é tomado como nulo (r = 0). Figura 6.2- Sistema planta mais controlador. 97 Representado o sistema da Figura (6.2) na forma de variáveis de estado, temos ẋ = Ax + Bu + Γw (6.20) y = Cx + v (6.21) onde x é o vetor de estado, u é o vetor de controle e y é o vetor de saı́das corrompidas por v; w e v são modelados como ruı́dos brancos, caracterizando processos estocásticos gaussianos com média zero. Considera-se que w e v não são correlacionadas no tempo e entre si, tendo as covariâncias: E wwT = Rf ≥ 0 E νν T = Qf > 0 E wν T = 0 (6.22) No problema LQG deseja-se minimizar a função custo: Z Jp = T T x Qc x + uT Rc u (6.23) 0 Onde as matrizes Qc e Rc são as mesmas definidas em 6.2.1. Considera-se que o vetor de saı́da seja os estados do sistema, nesse caso a matriz de saı́da é uma matriz identidade, ou seja: y = Cx (6.24) onde C = I, I é a matriz indentidade de ordem 6. O sistema 6.21 está sujeito a perturbações da planta e ruı́dos na observação da saı́da, a filosofia do projeto do controlador Kc pode ser estruturada (Athans, 1971) em três passos: 98 • projeto/análise determinı́stica do problema do controle; • projeto/análise do problema de estimação estocástica do estado; • projeto de um sistema de controle estocástico. A Figura (6.3) mostra a configuração do sistema de controle utilizado, note que o controle LQG é baseado em observador. Figura 6.3- Estrutura básica do sistema de controle LQG. A solução do problema LQG é conseguida pelo princı́pio da separação que possibilita a separação do problema original em dois problemas: Problema 1: Obter n o uma estimativa ótima x̂ do estado x de modo que T E (x̂ − x) (x̂ − x) seja minimizado. Esse sub-problema é resolvido pelo uso de um filtro de Kalman-Bucy, ignorando-se completamente o problema de controle. Problema 2: Obter o controlador para o problema linear quadrático determinı́stico (LQR considerando a referência zero (r = 0)), fazendo uso da estimativa x̂ como se fosse a medida exata do estado, ignorando-se completamente os aspectos estocásticos do problema. O problema do LQR é resolvido na seção 6.2.1. Nota: para a existência e estabilidade da solução do problema de estimação e/ou observação, a condição necessária e suficiente é que o sistema seja completamente 99 observável (Dorato e Cerone, 1995 e Maciejowski, 1989). As análises de observabilidade, mostradas em Kwakernaak e Sivan (1972), feitas para os sistema envolvido, Equação (5.33), garantem a solução do problema 2. A solução do problema 1 consiste na determinação da matriz ganho Kf . A estrutura do filtro de Kalman-Bucy é mostrado na Figura (6.4). Figura 6.4- Estrutura do filtro de Kalman-Bucy. A matriz ganho de Kalman Kf é dada por (Athans, 1971 e Maciejowski, 1989) Kf = Pf CT R−1 f (6.25) Onde Pf é dada pela Equação matricial algébrica de Ricatti 0 = −Pf A − AT Pf + Pf CT R−1 f Pf − Qf (6.26) É possı́vel mostrar que os autovalores do sistema completo, filtro mais controlador, são compostos pelo soma dos autovalores do filtro e do LQR (Kwakernaak e Sivan, 1972). As referências Flora (), Dorato e Cerone (1995) e Maciejowski (1989) mos100 tram que os autovalores do filtro e do LQR estão no semi-plano esquerdo (SPE) do plano complexo, ou seja os dois projetos são assintoticamente estáveis, então o sistema completo é, também assintoticamente estável. Na análise de robustez os dois projetos quando considerados isoladamente apresentam boas caracterı́sticas de robustez. Entretanto, como mostrado por Doyle (1978), o sistema completo, formado pelo filtro de Kalman-Bucy e pelo controlador determinı́stico não mantém as boas caracterı́sticas de robustez apresentadas, isoladamente, pelo filtro e pelo controlador. A recuperação das propriedades de robustez é conseguida pelo método LQG/LTR (loop transfer recovery), esse método é discutido em Maciejowski (1989) e Kwakernaak e Sivan (1972). 6.2.3 Controladores Baseados em Energia Nesta seção os controladores com realimentação de velocidade angular e atitude são investigados para o modo de estabilização. Ambos são baseados em energia e usam os atuadores magnéticos para o controle. Diz-se que esses controladores são baseados em energia pois a partir do cálculo da energia do sistema é escolhido uma lei de controle que satisfaz os critérios de estabilidade, usando para isso a teoria de Lyapunov. A função de Lyapunov V é escolhida como sendo a energia do sistema. As discussões são baseadas nas referencias Wisniewski (1996) e Fauske (2002). A prova da estabilidade dos controladores propostos e a obtenção dos pontos de equilı́brio estáveis do sistema são obtidos a partir da teoria de Lyapunov, e não fazem parte da proposta desse trabalho. Essas provas são fornecidas em Wisniewski (1996), Fauske (2002) e Overby (2004). Realimentação de Velocidade Angular A lei de controle com realimentação de velocidade angular é dada por (Fauske, 2002) mb = hω bob × Bb (6.27) onde h é uma constante escalar maior que zero (h > 0), que faz com que o sistema satélite equipado com bobinas, seja assintoticamente estável sobre quatro pontos de equilı́brio. 101 ω bob , cb3 , cb2 : (0, ±co3 , ±co2 ) (6.28) onde coi i = 1, 2, 3 são matrizes coluna tendo como componentes os cosenos diretores dos eixos do referencial do corpo em relação aos eixos do referencial orbital. Os quatro pontos de equilı́brio correspondem a orientação dos eixos x, y e z na mesma direção dos eixos xo , yo e zo , respectivamente. Realimentação de Atitude A lei de controle com realimentação de atitude é dada por (Fauske, 2002) mb = hω bob × Bb − α × Bb (6.29) onde é a componente vetorial do quaternion, Equação (4.7), e α é uma constante escalar, essa lei de controle faz com que o sistema, satélite equipado com bobinas, seja assintoticamente estável sobre a referência ou a VLHL. 6.3 ω bob , cb3 , cb2 : (0, co3 , co2 ) (6.30) Modo de Aquisição Para o estudo de aquisição de uma atitude arbitrária, definido nesse trabalho como modo de aquisição é feita a implementação da metodologia LQR rastreio (tracking) e um controlador proporcional derivativo (PD). O objetivo é avaliar, através da formulação desenvolvida nesse trabalho, a factibilidade do uso das rodas da TELDIX especificadas para o satélite EQUARS, para manobras de atitude e também comparar e discutir as leis de controle formuladas. Nas simulações os controladores são aplicados à dinâmica completa, não linear. 102 6.3.1 LQR Tracking Na aquisição de uma atitude arbitrária, onde a referência difere de zero (x 6= 0), o problema de formular uma lei de controle utilizando a metodologia LQR é conhecido como LQR rastreio/tracking. No LQR traking a lei de controle envolve um termo antecipativo (feedfoward ) adicionado ao termo de realimentação do estado (state-feedback ) (Dorato e Cerone, 1995 e Kirk, 1970), apresentado na seção 6.2.1. O problema de LQR tracking consiste em minimizar um ı́ndice de desempenho definido como Z Jp = T T x̃ Qc x̃ + uT Rc u (6.31) 0 onde x̃(t) = x(t) − xd (t) (6.32) onde xd é a trajetória de referência ou o estado desejado, x é obtido através do sistema dinâmico ẋ = A(t)x + B(t)u (6.33) No problema assumimos que o estado desejado (xd ) é conhecido e que o todos os estados x são disponı́veis para realimentação. As matrizes Qc e Rc são as mesmas definidas na seção 6.2.1. A lei de controle ótima para o problema de tracking é dada por (Kirk, 1970) u = −Kc (t)x + fw (t) (6.34) onde Kc (t) é a matriz de ganhos para realimentação dos estados, obtido pela Equação (6.8), obtido através da Equação de Ricatti, Equação (6.6), definida na seção 6.2.1: 103 Ṗ = −PA(t) − A(t)T P + PB(t)R−1 c P − Qc (6.35) e o termo antecipativo (fw (t)) é denominado sinal de comando e é dado por fw (t) = −Rc BT s(t) (6.36) onde s(t) é cálculado a partir da Equação diferencial linear (Kirk, 1970) −1 s(t) + Qc xd ṡ = − A(t) − PB(t)R−1 c B(t) (6.37) Note que o sinal de comando (fw (t)) depende dos parâmetros do sistema e do sinal de referência. A Figura (6.5) mostra a configuração do LQR tracking. Figura 6.5- Configuração do controle LQR tracking. 6.3.2 Proporcional Derivativo O controlador proporcional derivativo (PD), também é projetado a partir da planta linearizada, Equação (5.33). Nesse trabalho o projeto do controlador PD é baseado em Kaplan (1976). Uma discussão pormenorizada dos controladores industriais PD e PID é encontrada em Ogata (1998). Primeiro é obtida uma lei de controle para o movimento em torno do eixo de pitch, como se fosse um sistema 104 SISO, pois o movimento em pitch (θ) é desacoplado dos outros movimentos, roll e yaw. Em seguida é projetada uma lei de controle para o sistema de equações diferenciais acopladas, que descrevem o movimento de roll e yaw, e constitue um sistema MIMO com duas entradas e duas saı́das. Projeto em Pitch A partir da Equação (5.33), obtemos a Equação linearizada do movimento em pitch dada por J2 θ̈ + 3ωo2 (J1 − J3 ) θ + u2 = 0 (6.38) onde u2 representa o torque de controle da roda de reação disposta na direção do eixo de pitch. Note que a Equação (6.38) indica que o sistema não apresenta amortecimento, isso pode ser obtido através do sinal de controle u2 . Uma forma satisfatória para a lei de controle, discutida em Kaplan (1976) é dada por u2 = Kp2 Td2 θ̇ + θ − θr (6.39) onde o termo θ̇ introduz amortecimento ao sistema, Kp2 é o ganho proporcional (P) e Td2 a constante de tempo. O ganho derivativo é dado por Kd2 = Kp2 Td2 , θr é o ângulo de pitch de referência, considerado constante. Com a lei de controle o movimento em pitch torna-se o clássico sistema amortecido de segunda ordem. J2 θ̈ + Kp2 Td2 θ̇ + Kp2 + 3ωo2 (J1 − J3 ) θ − Kp2 θr = 0 (6.40) Os ganhos do controlador PD (Kp2 e Kd2 ), são obtidos nesse trabalho a partir da metodologia LQR, ou seja, a partir da matriz Kc tomando-se os ganhos da linha 2, que correspondem ao sinal de controle u2 . Uma outra alternativa para escolha dos ganhos é ajustar/sintonizar os ganhos a partir do diagrama de lugar das raı́zes, esse procedimento é usado em Kaplan (1976). 105 Projeto em Roll e Yaw A partir da Equação (5.33), obtemos as equações linearizadas do movimento em roll e yaw dadas por J1 φ̈ + a1 φ + a2 ψ̇ + u1 = 0 (6.41) J3 ψ̈ + c1 ψ + c2 φ̇ + u3 = 0 (6.42) onde os coeficientes a1 , a2 , c1 e c2 , são a1 = 4ωo2 (J2 − J3 ) + ωHo (6.43) a2 = −ωo (J1 − J2 + J3 ) + Ho (6.44) c1 = ωo2 (J2 − J1 ) + ωHo (6.45) c2 = ωo (J1 − J2 + J3 ) − Ho (6.46) os sinais de controle u1 e u3 são os torques de controle gerados pelas rodas de reação dispostas na direção do eixo de roll e yaw, respectivamente. Devido ao acoplamento dos movimentos as leis de controle são escolhidas como sendo u1 = Kp1 Td1 φ̇ + Kp1 (φ − φr ) + Kc1 ψ̇ (6.47) é lei de controle para o movimento em roll. O termo φ̇ introduz amortecimento ao sistema, Kp1 é o ganho proporcional (P), Td1 a constante de tempo. O ganho derivativo é dado por Kd1 = Kp1 Td1 , φr é o ângulo de roll de referência, considerado constante. O termo Kc1 é inserido devido ao acoplamento entre os movimentos em roll e yaw. Com a lei de controle o movimento em roll torna-se um sistema amortecido de segunda ordem. 106 J1 φ̈ + a1 φ + a2 ψ̇ + Kp1 Td1 φ̇ + Kp1 (φ − φr ) + Kc1 ψ̇ = 0 (6.48) Anologamente, a lei de controle para o movimento em yaw é escolhida como u3 = Kp3 Td3 ψ̇ + Kp3 (ψ − ψr ) + Kc3 φ̇ (6.49) O termo ψ̇ introduz amortecimento ao sistema, Kp3 é o ganho proporcional (P), Td3 a constante de tempo, o ganho derivativo é dado por Kd3 = Kp3 Td3 , ψr é o ângulo de yaw de referência, considerado constante. O termo Kc3 é inserido devido ao acoplamento entre os movimentos. Com a lei de controle o movimento em yaw torna-se em um sistema amortecido de segunda ordem. J3 ψ̈ + c1 ψ + c2 φ̇ + Kp3 Td3 ψ̇ + Kp3 (ψ − ψr ) + Kc3 φ̇ = 0 (6.50) Reescrevendo a equações (6.41) e (6.42) temos J1 φ̈ + Kp1 Td1 φ̇ + (a1 + Kp1 ) φ − Kp1 φr + (Kc1 + a2 ) ψ̇ = 0 (6.51) J3 ψ̈ + Kp3 Td3 ψ̇ + (c1 + Kp3 ) ψ − Kp3 ψr + (Kc3 + c2 ) φ̇ = 0 (6.52) A Figura (6.6) mostra a configuração do controlador PD. Os ganhos são obtidos a partir da metodogia LQR (matriz Kc ), assim como feito para o projeto de controle em pitch. A partir da matriz Kc toma-se os ganhos das linhas 1 e 3, que correspondem ao sinal de controle para u1 e u3 , respectivamente. Os valores de referência (φr , θr , ψr ) correspondem a uma atitude nominal arbitrária. A vantagem do uso do controlador PD é que a realimentação é apenas dos ângulos de atitude, não sendo necessário medir as taxas de variação dos ângulos, entretanto, o procedimento de derivação dos ângulos gera erros. 107 Figura 6.6- Configuração do controlador PD. 108 CAPÍTULO 7 SIMULAÇÕES Neste Capı́tulo, são feitas as simulações utilizando as estratégias de controle para os modos de operação: • Detumble; • Estabilização; • Aquisição. As configurações propostas para o ACS utilizam como atuadores: rodas de reação combinadas com bobinas magnéticas e apenas bobinas magnéticas para o controle autônomo do satélite, nos dois primeiros modos de operação (detumble e estabilização). O controlador Bdot é utilizado para o descapotamento/detumbling do satélite. O controlador proposto para estabilização e manutenção da atitude nominal é obtido usando a metodologia LQR e LQG. Outra alternativa apresentada para fase de estabilização e manutenção da atitude nominal, usando bobinas magnéticas, são os controladores baseados em energia: realimentação da velocidade angular (angular velocity feedback ) e realimentação da atitude ()(attitude feedback ). Nas simulações são usados os parâmetros do satélite brasileiro EQUARS. Para a fase de aquisição é analisada a exeqüibilidade da utilização das rodas de reação onde é proposto o controle PD e LQR rastreio/tracking. Os parâmetros usados na simulação são mostrados na Tabela (7.1). Os parâmetros e especificações obtidos para o satélite EQUARS são baseados em Carvalho (2003) e Heidelberg (2004). As três rodas de reação são iguais e as três bobinas também são iguais. 7.1 Modo de Detumble O controlador usado para o modo de detumble foi o controlador Bdot. As condições iniciais e os parâmetros do controlador são listados a seguir. 109 Tabela 7.1- Parâmetros de simulação. Parâmetros Valores Jx = 13.31, Jy = 14.22, Jz = 11.20 kgm2 1o sobre cada eixo (φ, θ, ψ) Tensor de inércia do satélite (J) Precisão de apontamento Bobinas magnéticas: Máximo momento dipolo Seção da área Número de espiras Resistência 2 Am2 0.075 m2 100 20 Ω Rodas de reação: Momento de inércia Máxima rotação por minuto Máxima quantidade de movimento angular Máximo torque Parâmetros orbitais: Altura (h) Excentricidade (e) Ascenção do nodo ascendente (Ω) Inclinação (i) Velocidade orbital média Perı́odo orbital Data da simulação Jw = 0.015 kgm2 7500 rpm 12 N ms 75 mN m 750 km ∼ =0 30 graus 20 graus ωo = 0.001 rad/s To ∼ = 1h40mim 31 − 12 − 2004 110 Condições iniciais Parâmetros de controlador b ωib = [10 10 10]T graus/s k = 2 · 105 , mc = −0.1 Am2 onde a lei de controle Bdot é dada pela Equação (6.1) mb = −k Ḃb − mc (7.1) A escolha da velocidade angular inicial foi baseada nas simulações feitas para o modo de estabilização, Seção (6.2). Para essa condição inicial de velocidade angular b b b (ωix , ωiy , ωiz ) (≈ 2 rpm) os atuadores especificados saturam rapidamente no modo de estabilização e não são capazes de estabilizar ou controlar a atitude do veı́culo, definindo uma situação ou velocidade angular crı́tica. Os parâmetros do controlador foram obtidos por tentativa. A Figura (7.1) mostra a redução da velocidade angular do satélite. Depois de quatro revoluções em volta da Terra, o veı́culo já está apto a executar manobras de atitude, ou seja, entrar no modo de estabilização. Figura 7.1- Velocidade angular do satélite ω bib . A Figura (7.2) mostra o amortecimento da velocidade angular do satélitea partir da quarta órbita. Reduzindo a velocidade de rotação do satélite para abaixo de (2 graus/s) . 111 Figura 7.2- velocidade angular ω bib . A Figura (7.3) mostra o dipolo magnético das bobinas, o torque de controle gerado e a potência de cada bobina. Os valores do dipolo magnético atendem as especificações dos atuadores magnéticos (mx , my e mz < 2 Am2 ). Figura 7.3- Dipolo magnético, torque magnético e potência. 112 A Figura (7.4) mostra os mesmos resultados mas tomados a partir da quarta órbita. Figura 7.4- Dipolo magnético, torque magnético e potência. A lei de controle apresentada (Bdot) é de simples implentação e não requer informações da atitude do veı́culo. Apesar das simulações terem avaliado uma condição inicial crı́tica de velocidade angular, essa pode ser maior, dependendo de como o foguete lança o satélite na sua órbita, ou seja, depende do veı́culo lancador. Ainda assim o controlador Bdot pode ser utilizado, entretanto, exigindo um tempo maior para a redução da velocidade angular. 7.2 Modo de Estabilização As estratégias de controle usada no modo de estabilização são baseadas nas metodologias LQR e LQG, apresentadas no Capı́tulo 6. A metodologia LQR é implementada nos dois modelos: satélite equipado com bobinas e satélite equipado com rodas. A metodologia LQG é usada no modelo do satélite equipado com rodas. Os controladores baseados em energia são usados no modelo do satélite equipado com bobinas. 113 7.2.1 Método LQR Nessa Seção o método LQR é usado na formulação da lei de controle do modelo do satélite equipado com rodas. São investigados três cenários: • condições iniciais crı́ticas (ωx , ωy , ωz > 1 rpm) ; • condições iniciais quase crı́ticas (ωx , ωy , ωz = 1 rpm); • velocidade angular (ωx , ωy , ωz ∼ = 0 rpm). O objetivo é avaliar situações nas quais o projeto LQR ainda é válido no controle de atitude em três eixos. A Tabela (7.2) mostra as condições iniciais para os três casos. Tabela 7.2- Condições iniciais da simulação para o modo de estabilização utilizando a metodologia LQR. casos 1 2 3 velocidade angular inicial atitude inicial ω bib = [0 0 0]T graus/s ω bib = [6 6 6]T graus/s ω bib = [10 10 10]T graus/s (φ, θ, ψ) = (30, −20, 25) graus (φ, θ, ψ) = (30, −20, 25) graus (φ, θ, ψ) = (30, −20, 25) graus Em seguida será usada a formulação da lei de controle do modelo do satélite equipado com bobinas. O objetivo e estudar/analisar os dois procedimentos para o modo de estabilização Satélite Equipado com Rodas Caso 1: Velocidade Angular Para o primeiro caso considera-se que a velocidade angular do satélite seja nula, conseguida através do controlador Bdot. A partir dessa condição inicial temos que a lei de controle para o modo de estabilização é dada pela Equação (6.7), reescrita aqui 114 u = −Kc x (7.2) Kc = R−1 c BPc (7.3) onde onde Pc é dado pela solução da Equação de Ricatti, Equação (6.6). Os valores númericos das matrizes A e B, usadas no projeto de controle, são obtidos a partir das equações (5.34) e (5.35). As matrizes de ponderação Rc e Qc selecionadas são mostradas na Seção (6.2.1). Resolvendo a Equação de Ricatti, obtemos o ganho usado na lei de controle −0.0286 0.0000 0.0004 −0.8476 0.0000 0.0000 Kc = 0.0000 −0.0287 0.0000 0.000 −0.9180 0.0000 −0.0004 0.0000 −0.0286 0.000 0.0000 0.9510 (7.4) É importante resaltar que em todas as simulações aplicou-se a lei de controle linear para o modelo completo, não linear. A Figura (7.5), mostra os ângulos de atitude em função do tempo, durante o modo de estabilização. O resultado mostra um bom desempenho do sistema de controle. Note que a partir de 100 segundos o satélite já atinge a precisão requerida (< 1o ) para a atitude nominal especificada, (VLHL), realizando em seguida a manutenção da atitude nominal. A Figura (7.6) mostra o torque gerado pelas três rodas de reação. O primeiro gráfico mostra o torque (τ bw ) gerado pelas rodas dispostas nas três direções principais de inércia. O segundo gráfico mostra o torque de acoplamento devido ao movimento de rotação do satélite, sendo da ordem de 100 vezes menor que o torque τ bw gerado pelas rodas. O torque efetivo de controle é a soma dos dois torques. Para efeito de projeto da roda o efeito de acoplamento é considerado uma perturbação. O último gráfico mostra a quantidade de movimento angular das rodas, note que a roda ao longo do eixo de pitch (normal a órbita) tem uma quantidade de movimento angular nominal diferente de zero, fornecendo rı́gidez giroscópica ao sistema durante a manutenção 115 Figura 7.5- Ângulos de atitude roll, pitch e yaw em função do tempo para o caso 1. da atitude nominal. Figura 7.6- Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas para o caso 1. A Figura (7.7) mostra as rotações por minuto de cada uma das rodas e as velocidades 116 angulares ω bib e ω bob , respectivamente. Como era de se esperar a velocidade relativa entre os referenciais do satélite (BF) e orbital (OF) é nula. Figura 7.7- Velocidades angulares ω bib e ω bob e rotações por minuto das rodas de reação. A simulações mostram que, para a condição inicial de velocidade angular assumida, Tabela (7.2), todos os requisitos de precisão e operação são atendidos, não havendo saturação dos atuadores no modo de estabilização. Caso 2: Velocidade Angular Quase Crı́tica Para o segundo caso considera-se que a velocidade angular residual do satélite seja de 1 rpm, que equivale a 6 graus/s. A partir dessa condição inicial usa-se a lei de controle dada pela Equação (6.7), como no caso 1. O objetivo é avaliar se mesmo sob uma condição de velocidade angular não nula. Esta configuração poderia, por exemplo, ser causada por imprecisão de medidas ou falha no sistema de controle para o modo de detumble. Os resultados mostram que ainda assim, o sistema de controle ainda é capaz de estabilizar e controlar o movimento de atitude, a partir dos atuadores (rodas) especı́ficadas para o satélite EQUARS. Os valores de Kc obtidos e as matrizes de ponderação Rc e Qc selecionadas são as mesmas utilizadas no caso 1. 117 A Figura (7.8), mostra os ângulos de atitude durante o modo de estabilização. O resultado mostra um desempenho pobre do sistema de controle. Note que, apesar do desempenho, a partir de 200 segundos o satélite atinge a precisão requerida (< 1o ) para a atitude nominal especificada (VLHL), realizando em seguida manutenção da atitude nominal sem maiores problemas. Figura 7.8- Ângulos de atitude roll, pitch e yaw em função do tempo para o caso 2. A Figura (7.9) mostra o torque gerado pelas três rodas de reação, o primeiro gráfico mostra o torque (τ bw ) gerado pelas rodas, o gráfico mostra, claramente, a saturação dos atuadores, todas as rodas em algum momento sofrem saturação o segundo gráfico mostra o torque de acoplamento devido ao movimento de rotação do satélite, sendo bem maior a influência que no caso 1. O último gráfico mostra a quantidade de movimento angular das rodas, note que apesar da saturação em torques (> 75mN m), nenhuma das rodas apresentam saturação na quantidade de movimento angular (< 12N ms). A Figura (7.10) mostra as rotações por minuto de cada uma das rodas e as velocidades angulares ω bib e ω bob . A simulações mostram que para a condição inicial de velocidade angular assumida (Tabela (7.2)) mesmo sob saturação das rodas, todos os requisitos de precisão são 118 Figura 7.9- Torque τ bw , torque de acoplamento e quantidade de movimento angular das rodas para o caso 2. Figura 7.10- Velocidades angulares ω bib e ω bob e rotações por minuto das rodas para o caso 2 . atendidos pelo sistema de controle. 119 Caso 3: Velocidade Angular Crı́tica Para o terceiro caso analisado usamos o mesmo projeto de controle dos casos 1 e 2, e as condições iniciais mostradas na Tabela (7.2). A Figura (7.11) mostra os ângulos de atitude em relação ao tempo. Nota-se que o sistema de controle não é capaz de atender os requisitos no modo de estabilização, devido a uma velocidade angular crı́tica de aproximamente 1.5 rotações por minuto (rpm). Portanto, para o projeto do sistema de controle avaliado é necessário garantir que no modo de estabilização a velocidade angular inicial seja igual ou menor que 1.5 rpm. Essa redução é obtida no modo de detumble como mostrado nas simulações, Seção (7.1). Figura 7.11- Ângulos de atitude roll, pitch e yaw e tempo para o caso 3. A Figura (7.12) mostra os gráficos das velocidades angulares ω bib e ω bob e as rotações por minuto das rodas, nota-se que as rodas dispostas ao longo dos eixos de roll e yaw saturam (> 7500rpm). A roda em roll atinge rapidamente sua velocidade máxima de operação. A Figura (7.13) mostra a saturação em torque e quantidade de movimento angular das rodas. Os resultados da simulação mostram que o sistema de controle não é capaz de 120 Figura 7.12- Velocidades angulares ω bib e ω bob e rotações por minuto das rodas para o caso 3. Figura 7.13- Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas para o caso 3. estabilizar e controlar a atitude. Isso se deve, fundamentalmente, a limitação dos atuadores e a dinâmica se tornar altamente não linear. 121 Satélite Equipado com Bobinas A formulação da lei de controle utilizando a metodologia LQR para o modelo do satélite equipado com bobinas foi discutida na Seção (6.2.1). O ganho de controle Kc é variante no tempo, dada pela Equação (6.13) u = −Kc (t)x (7.5) A matriz solução da Equação de Ricatti Pc , Equação (6.6), é calculada para cada passo de integração, sendo aproximadamente periódica devido ao campo geomagnético (Bo ) ser quase periódico. Os valores númericos das matrizes A, B, usadas no projeto de controle, são obtidos a partir das equações (5.48) e (5.49). A Figura (7.14) mostra a variação do elemento Kc (1, 1) da matriz de ganho do controlador, mostrando claramente que os ganhos são aproximadamente periódicos. Isso permite reduzir o esforço computacional introduzindo uma matriz de ganhos previamente calculada para uma órbita nominal. Figura 7.14- Ganho do controlador em função das órbitas Kc (1, 1). As matrizes de ponderação Rc e Qc selecionadas são mostradas na Seção (6.2.1). As condições iniciais são listadas a seguir 122 Condições iniciais: Velocidade angular Atitude inicial ω bib = [0 0 0]T graus/s (φ, θ, ψ) = (15, −10, 10) Novamente utilizou-se nas simulações o modelo completo, não linear e no caso do satélite equipado com bobinas o modelo também é variante no tempo. A Figura (7.15) mostra os ângulos de atitude durante o modo de estabilização, o resultado mostra um desempenho razoável do sistema de controle, a partir dos parâmetros selecionados para o controlador. Note que é necessário duas revoluções do satélite em torno da Terra (≈ 3h20min) para que esse atinga a precisão requerida (< 1o ) realizando em seguida a manutenção da atitude nominal, ou seja, alinhado com a vertical local (VLHL). Figura 7.15- Ângulos de atitude roll, pitch e yaw em função do número de órbitas, utilizando bobinas. A Figura (7.16) mostra o dipolo magnético das bobinas, o torque de controle gerado, e a potência de cada bobina, os valores do dipolo magnético atendem as especificações dos atuadores magnéticos (mx , my e mz < 2 Am2 ), ficando bem abaixo do dimensionado para o sistema de controle de atitude. Apesar das bobinas apresentarem um tempo de resposta bem maior, a potência exigida pelo sistema de controle é muito menor do que aquela usada no sistema de 123 Figura 7.16- Dipolo magnético, torque de controle (τ bm ) e potência das bobinas em função do número de órbitas, utilizando bobinas. controle de atitude do satélite equipado com rodas; tornando-se uma opção atrativa, devido a seu custo e baixo consumo de energia. 7.2.2 Controladores Baseados em Energia Uma alternativa para o sistema de controle de atitude empregando bobinas magnéticas são os controladores não lineares propostos na Seção (6.2.3). Sua principal vantagem é a de que as leis de controle são formuladas a partir do sistema não linear, baseado nos crı́terios de estabilidade de Lyapunov. Isso garante a estabilização, não sendo válida somente em torno de um ponto de operação, como no caso de sistema lineares. São analisados duas leis de controle: realimentação de velocidade angular e realimentação de atitude. Realimentação de Velocidade Angular A lei de controle é dada pela Equação (6.27) 124 mb = hω bob × Bb (7.6) As condições iniciais e os parâmetros do controlador obtidos por tentativa são listados a seguir Condições iniciais: Velocidade angular ω bib = [0 0 0]T graus/s Atitude inicial (φ, θ, ψ) = (60, −30, 40)graus Parâmetros do controlador h = 3.25 · 106 A Figura (7.17) mostra os ângulos de atitude no modo de estabilização usando a realimentação de velocidade angular. O resultado mostra um desempenho razoável para os ângulos de roll e pitch. Para o eixo de yaw é obtido um pobre desempenho do sistema de controle. Esse resultado ocorre devido a falta de informação da atitude. É necessário quatro revoluções do satélite em torno da Terra (≈ 6h40min) para que esse atinga a precisão requerida (< 1o ) nos eixos de roll e pitch, sendo necessário quase o dobro de revoluções para o eixo de yaw. Em seguida realiza-se a manutenção da atitude nominal, ou seja, alinhado com a vertical local (VLHL). Figura 7.17- Ângulos de atitude em função do número de órbitas. A Figura (7.18) mostra o dipolo magnético das bobinas, o torque de controle gerado, 125 e a potência de cada bobina, os valores do dipolo magnético atendem as especificações dos atuadores magnéticos, ficando abaixo do dimensionado para o sistema de controle de atitude. O terceiro gráfico mostra que a potência exigida pelas bobinas é bem baixa. Figura 7.18- Dipolo magnético, torque de controle e potência das bobinas em função do número de órbitas. A Figura (7.19) mostra as velocidades angulares e o campo geomagnético expresso no sistema do satélite. É possı́vel mostrar (Fauske, 2002 e Overby, 2004) que a lei de controle apresentada não garante, sobre outras condições iniciais, a estabilização e o controle em três eixos de interesse, ou seja, alinhado com a VLHL. Isso se deve, ao sistema dinâmico ter mais de um ponto de equilı́brio, sendo necessário informações da atitude para a estabilização de interesse. A lei de controle apresentada é uma alternativa que deve ser considerada, para os casos em que nenhuma informação da atitude é disponı́vel para o sistema de controle, devido a possı́veis falhas nos sensores de atitude. Entretanto, não se pode garantir, em alguns casos, o controle em para a orientação de interesse (VLHL). Novamente o uso de bobinas apresenta um baixo consumo de energia. 126 Figura 7.19- Velocidades angulares ω bib e ω bob e campo magnético terrestre local Bb . Realimentação de Atitude A lei de controle é dada pela Equação (6.29) mb = hω bob × Bb − α × Bb (7.7) As condições iniciais e os parâmetros do controlador são listados a seguir Condições iniciais: Velocidade angular ω bib = [0 0 0]T graus/s Atitude inicial (φ, θ, ψ) = (60, −30, 40)graus Parâmetros do controlador h = 3.25 · 106 α = −3000 A Figura (7.20) mostra os ângulos de atitude no modo de estabilização usando a realimentação de velocidade angular e atitude. O resultado mostra um desempenho razoável para os ângulos de roll, pitch e yaw. É necessário duas revoluções do satélite em torno da Terra (≈ 3h20min) para que esse atinja a precisão requerida (< 1o ) nos eixos de roll, pitch e yaw. O desempenho do sistema de controle é próximo 127 ao obtido com a metodologia LQR, entretanto, a atitude inicial assumida é bem maior. O resultado mostra que a informação da atitude é importante para um bom desempenho do sistema de controle, na estabilização da atitude em três eixos. A partir da segunda revolução do satélite em torno da Terra ocorre a manutenção da atitude nominal (VLHL), ou seja, o sistema se estabiliza em torno da atitude nominal. Figura 7.20- Ângulos de atitude. A Figura (7.21) mostra o dipolo magnético das bobinas, o torque de controle gerado, e a potência de cada bobina. Os valores do dipolo magnético atendem as especificações dos atuadores magnéticos, Tabela (7.1), ficando abaixo do dimensionado para o sistema de controle de atitude. O terceiro gráfico mostra que a potência exigida pelas bobinas é bem baixa. A Figura (7.22) mostra as velocidades angulares e o campo geomagnético expresso no sistema do satélite. A lei de controle apresentada é uma boa alternativa para o sistema de controle de atitude em três eixos, pois apresenta um desempenho razoável. A principal vantagem em relação a metodologia LQR é que mesmo sob as não linearidades, a eficiência do sistema de controle de atitude na estabilização do sistema é garantida. Novamente o 128 Figura 7.21- Dipolo magnético, torque de controle e potência das bobinas em função do número de órbitas. Figura 7.22- Velocidades angulares ω bib e ω bob e campo magnético local Bb em função do número de órbitas. uso de bobinas apresenta uma opção atrativa devido a baixa potência, o que equivale a um baixo consumo de energia. 129 7.2.3 Método LQG Como descrito na Seção (6.2.2), metodogia LQG é usada no modo de estabilização para o satélite equipado com rodas. As condições iniciais, e os parâmetros do filtro de Kalman são listados na Tabela (7.3) Tabela 7.3- Parâmetros de simulação. Condições iniciais: Velocidade angular ω bib Atitude inicial (φ, θ, ψ) Parâmetros do filtro: Matriz Rf Matriz Qf [0 0 0]T graus/s (30, −20, 25)graus diag[(1, 1, 1, 0.1, 0.1, 0.1)]0.52 × (π/180o ) diag[(1, 1, 1, 1, 1, 1)]152 × (π/180o ) A matriz Rf é obtida a partir do desvio padrão de sensores solares (0.5o ) e de girômetros (0.05graus/s) tı́picos. A matriz Q é obtida através de ajustes por tentativa. A partir dos parâmetros do filtro obtemos o matriz de ganho de Kalman 30.0002 0.0000 0.0000 0.9091 0.0000 0.0000 0.0000 30.0002 0.0000 0.000 0.9109 0.0000 0.0000 0.0000 30.0002 0.000 0.0000 0.9509 Kf = 0.0000 0.0000 0.0000 300.0000 0.0000 0.0001 0.0000 0.0091 0.0000 0.0000 300.0000 0.0000 0.0000 0.0000 0.0091 0.001 0.0000 300.0000 (7.8) Os parâmetros do controlador usado são os mesmos utilizados no método LQR. A Figura (7.23) mostra os ângulos de atitude ou de Euler estimados com o filtro de Kalman. O desempenho do controlador é muito próximo ao obtido com o projeto LQR (caso 1). A Figura (7.24), mostra os ângulos de Euler simulados, esses são corrompidas por ruı́dos de distribuição aleatória e uniforme. 130 Figura 7.23- Ângulos de atitude estimados em função do tempo. Figura 7.24- Ângulos de atitude simulados. A Figura (7.25) confronta os ângulos de atitude estimados e simulados, mostrando claramente a suavização do sinal estimado, que é usado na realimentação do sistema de controle. A Figura (7.26) mostra os gráficos do erro (valor real (simulado) − valor estimado). O primeiro gráfico mostra o erro nos ângulos de atitude e o segundo 131 gráfico mostra o erro nas taxas de variação dos ângulos de atitude. Em ambos os casos o erro é menor que a covariância em ∼ = 70% dos casos, representada pela linha vermelha. Figura 7.25- Ângulos de atitude simulados e estimados. Figura 7.26- Erro dos ângulos de atitude e variação da atitude. A Figura (7.27) mostra o torque gerado pelas três rodas de reação, o torque de 132 acoplamento e a quantidade de movimento angular. O primeiro gráfico mostra o torque (τ bw ) gerado por cada roda, o segundo gráfico mostra o torque de acoplamento devido ao movimento de rotação do satélite. As magnitudes ou esforço de controle são próximas às obtida no projeto LQR (caso 1), como era de se esperar, pois o controlador é o mesmo desenvolvido no caso 1. O último gráfico mostra a quantidade de movimento angular das rodas. A roda ao longo do eixo de pitch (normal a órbita) tem uma quantidade de movimento angular nominal diferente de zero, fornecendo rı́gidez giroscópica ao sistema durante a manutenção da atitude. Figura 7.27- Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas. 7.3 Modo de Aquisição Para o sistema de controle no modo de aquisição são avaliados duas alternativas para o projeto de controle: 1) controlador Proporcional Derivativo (PD); 2) LQR rastreio/tracking. Esse projeto tem o objetivo de avaliar a factibilidade do emprego das rodas especificadas para o satélite brasileiro EQUARS (ver Heidelberg, 2004), na execusão de manobras de grandes ângulos. No modo de aquisição, assumese que a condição inicial o satélite está alinhado com a vertical e horizontal local (VLHL), ou seja, o satélite parte do modo de estabilização. As condições iniciais e a atitude de refêrencia são mostradas na Tabela (7.4). 133 Tabela 7.4- Condições iniciais da simulação para o modo de aquisição utilizando rodas. Lei de controle velocidade angular inicial atitude inicial (φ, θ, ψ) atitude de referência (φr , θr , ψr ) LQR tracking PD ω bib = [0 0 0]T graus/s ω bib = [0 0 0]T graus/s (0, 0, 0) graus (0, 0, 0) graus (−60, 70, 130) graus (−60, 70, 130) graus 7.3.1 LQR Tracking A lei de controle ao método LQR tracking no modo de aquisição, apresentada na Seção (6.3.1), é dada pela Equação (6.34), u = −Kc (t)x + fw (t) (7.9) onde a matriz Kc é obtida através da solução da Equação de Ricatti. As matrizes de ponderação Qc e Rc são as mesma utilizadas no projeto LQR. O sinal de comando é dado pela Equação (6.36). fw (t) = −Rc BT s(t) (7.10) onde o vetor s(t) é função da referência, atitude (φr , θr , ψr ) e velocidades angulares (φ̇r , θ̇r , ψ̇r ), calculada pela Equação (6.37). As matrizes obtidas no projeto LQR tracking são: −0.0286 0.0000 0.0004 −0.8476 0.0000 0.0000 Kc = 0.0000 −0.0286 0.0000 0.0000 −0.9180 0.0000 0.0004 0.0000 −0.0286 0.0000 0.0000 −0.9510 e 134 (7.11) 0.1015 −0.1285 −0.2472 s= 1.3734 −1.8581 −3.7070 (7.12) A Figura (7.28) mostra os ângulos de atitude durante o modo de aquisição. O resultado mostra um bom desempenho do sistema do controle, usando o método LQR tracking. A partir de 200 segundos o satélite atingi a precisão requerida (< 1o ) para a atitude especificada ou de referência (φr , θr , ψr ) = (−60, 70, 130). O modelo do satélite usado nas simulações é o modelo completo, não linear. Figura 7.28- Ângulos de atitude roll, pitch e yaw em função do tempo para a aquisição de atitude. A Figura (7.29) mostra o torque gerado pelas três rodas de reação, o primeiro gráfico mostra o torque (τ bw ) gerado pelas rodas dispostas nas três direções principais de inércia o segundo gráfico mostra o torque de acoplamento devido ao movimento de rotação do satélite. O torque efetivo de controle é a soma dos dois torques. O resultado mostra que mesmo para manobras de grandes ângulos, não há saturação da quantidade de movimento das rodas (hw < 12N ms), ficando abaixo do valor 135 de saturação. O primeiro grafico mostra que não ocorre saturação do esforço ou torque de controle (τ bw < 75mN m). Esse resultado mostra a viabilidade da utilização das rodas especificadas para o satélite EQUARS para a aquisição de uma atitude arbitrária. Figura 7.29- Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas em função do tempo. A Figura (7.30) mostra as rotações em rpm de cada uma das rodas e as velocidades angulares ω bib e ω bob . Na manutenção da atitude de referência as velocidades das rodas são nominalmente diferentes de zero. Nas simulações, apesar do projeto do controlador LQR tracking se basear na planta linear, o controle empregado sobre a planta não linear apresenta boas caracterı́sticas de desempenho. O resultado se deve em parte as boas propriedades de robustez do controle LQR. O controle LQR é uma alternativa atraente devido a otimalidade. A escolha de ponderações para o estado e controle, permite ao projetista um projeto de controle com base em informações e especificações dos atuadores e respectivas faixas de operação. Um menor esforço para o ajuste dos ganhos Qc e Rc também é conseguido, através do algoritmo apresentado, Equação (6.17). 136 Figura 7.30- Rotações por minuto das rodas de reação e velocidades angulares ω bib e ω bob . 7.3.2 Controle PD Os ganhos para o controlador PD foram obtidos usando a metodologia LQ, Seção (6.2.1). A lei de controle PD, obtida na Seção (6.3.2), equações (6.47), (6.39) e (6.49) podem ser reescritas como φ − φr φ̇ u = Kp θ − θr + Kd θ̇ ψ̇ ψ − ψr (7.13) Onde as matrizes Kp e Kd são as matrizes com os ganhos proporcionais e derivativos, respectivamente. Os ganhos do controlador PD obtidos são −0.0286 0.0000 0.0000 Kp = 0.0000 −0.0286 0.0000 0.0000 0.0000 −0.0286 137 (7.14) e −0.8476 0.0000 6.2714 · 10−5 Kd = 0.0000 −0.9180 0.0000 8.1064 · 10−5 0.0000 −0.9510 (7.15) A Figura (7.31) mostra os ângulos de atitude durante o modo de aquisição. O resultado mostra um bom desempenho do sistema do controle, a partir dos ganhos Kp e Kd usados para o controlador PD. A partir de 200 segundos o satélite já atinge a precisão requerida (< 1o ) para a atitude especificada ou de referência (φr , θr , ψr ) = (−60, 70, 130). Novamente foi usado o modelo completo (não linear), nas simulações. Figura 7.31- Ângulos de Euler em função do tempo. A Figura (7.32) mostra o torque gerado pelas três rodas de reação. O primeiro gráfico mostra o torque (τ bw ) gerado pelas rodas dispostas nas três direções principais de inércia. O segundo gráfico mostra o torque de acoplamento devido ao movimento de rotação do satélite. O torque efetivo de controle é a soma dos dois torques. O último gráfico mostra a quantidade de movimento angular das rodas. Note que as três rodas ao longo dos eixos de roll, pitch e yaw apresentam uma quantidade de 138 movimento angular nominal diferente de zero, fornecendo rı́gidez giroscópica ao sistema durante a manutenção da atitude de referência. O resultado mostra que mesmo para manobras de grandes ângulos não há saturação de quantidade de movimento angular das rodas (hw < 12N ms), ficando abaixo desse valor. O primeiro grafico mostra que não ocorre saturação do esforço de controle (τ bw < 75mN m), chegando a um valor máximo no eixo de yaw de ∼ = 60mN m. Esse resultado mostra que as rodas especificadas para o satélite EQUARS podem ser usadas para realizar manobras de grandes ângulos, ou seja, ser usadas para aquisição de uma atitude arbitrária. Figura 7.32- Torque τ bw , torque de acoplamente e quantidade de movimento angular das rodas em função do tempo. A Figura (7.33) mostra as rotações por minuto de cada uma das rodas e as velocidades angulares ω bib e ω bob . Na manutenção da atitude de referência as velocidades das rodas são nominalmente diferentes de zero. O resultados obtidos mostram que é factı́vel o uso das rodas especificadas (Tabela 7.1), para o modo de aquisição do satélite EQUARS. Uma das vantagens no uso do controlador PD em relação ao LQR é usar como realimentação apenas a atitude. Entretanto, tem como incoveniente a derivação dos sinais. A metodologia LQ usada para seleção/ajuste dos ganhos mostra-se muito prática para o controle PD. 139 Figura 7.33- Rotações por minuto das rodas e velocidades angulares ω bib e ω bob em função do tempo. 140 CAPÍTULO 8 CONCLUSÃO Nesse Capı́tulo apresentam-se as principais conclusões relacionadas aos resultados obtidos, encerrando com sujestões para trabalhos futuros. Nesse trabalhos apresentou-se um estudo de alternativas de sistemas de controle (ACSs) para satélites artificiais estabilizados em três eixos. Os procedimentos adotados foram: 1) satélite utilizando rodas e bobinas como atuadores e 2) satélite utilizando apenas bobinas magnéticas. Para cada um dos procedimentos analisou-se a exeqüibilidade dos ACS(s) para os dois principais modos de operação estudados: 1) detumble, 2) estabilização. As leis de controle empregadas para esses modos foram LQR, LQG e controladores baseados em energia (Bdot, realimentação de atitude e velocidade angular). Os métodos LQR e LQG resultam em projetos lineares e os reguladores baseados em energia resultam em projetos não lineares. Apesar do projeto linear os métodos LQR e LQG foram aplicados à dinâmica não linear, apresentando bons resultados, no que se refere ao desempenho e robustez, também fornecendo boas margens de ganho e de fase. A Tabela (8.1) resume as diferentes configurações dos ACS(s) estudados e as diferentes estratégias de controle. Tabela 8.1- Alternativas para o ACS. Modo ACS Lei de controle Detumble bobinas Rodas e bobinas Bdot Bdot Estabilização Rodas e bobinas Rodas e bobinas Bobinas Bobinas Bobinas LQR LQG LQR Realimentação de velocidade angular Realimentação de atitude Aquisição Rodas e Bobinas PD Rodas e Bobinas LQR tracking A lei de controle Bdot utilizada no modo de detumble, apresentou um bom resul141 tado em termos de tempo, reduzindo a velocidade de rotação do satélite para um nı́vel tı́pico admissı́vel (< 0.01rpm), colocando o satélite em condições de inicializar o modo de estabilização. A partir de uma velocidade crı́tica ∼ = 2rpm(definida nesse trabalho como a velocidade em que o ACS no modo de estabilização não é apropriado, devido a saturação dos atuadores) a estratégia de controle leva aproximadamente quatro órbitas para atingir uma faixa de velocidades angulares admissı́veis. A estratégia de controle é factı́vel de ser empregada, tomando como base a consumo de energia dos atuadores e valores de saturação. A análise de factibilidade do sistema de controle foi feita com base na potência especificada para satélites tı́picos do EQUARS. A lei de controle é fácil de ser implementada e não requer informações ou conhecimento da atitude do veı́culo, mas apenas dados dos magnetômetros. Apesar das simulações terem avaliado uma condição inicial de velocidade angular crı́tica de ∼ = 2rpm, ela pode ser maior dependendo do veı́culo lançador. Entretanto, o controlador Bdot ainda pode ser usado, se o requisito de tempo não for muito estreito, pois exigiria um tempo maior para a redução da velocidade angular. Para o modo de estabilização foram comparados o uso de bobinas e rodas. Na utilização de rodas foram analisadas três cenários: 1) com velocidade residual de rotação nula (ω bib = [0, 0, 0]T graus/s), 2) com velocidade quase crı́tica (ω bib = [6, 6, 6]T graus/s) e 3) com velocidade crı́tica (ω bib = [10, 10, 10]T graus/s). Os resultados mostram que para o primeiro caso o sistema de controle consegue atender os requisitos de operação muito bem, não ocorrendo saturação na velocidade da rodas nem em esforço de controle ou torque. Para o segundo caso mesmo ocorrendo saturação no torque o sistema de controle consegue estabilizar a atitude e realizar a manutenção ou o alinhamento com a VLHL. Para o terceiro caso o sistema de controle desestabiliza a atitude devido a saturação nos atuadores. Uma solução para esse problema seria ajustar os ganhos do controlador através das matrizes de ponderação do estado e controle (Qc e Rc ), Diminuindo o esforço de controle e aumentado o tempo para estabilização que nos casos 1 e 2 estudados é relativamente rápido (< 200s). O método LQR desenvolvido para os dois modelos: satélite equipado com rodas e satélite equipado com bobinas, atende as especificações do sistema de controle em ambos os casos. As principais caracterı́sticas apresentadas pelo método são: controle ótimo, robustez, fácil implementação, ponderações no estado e controle permitindo ao projetista um maior sentimento fı́sico do problema. A principal desvantagem do 142 método é o esforço computacional para o caso de sistema variantes no tempo, como é o caso do satélite equipado com bobinas. Entretanto isso pode ser contornado estocando-se um ganho previamente calculado que, devido a quase periodicidade do campo geomagnético, o ganho é quase periódico, com mostrado nas simulações. Outro importante resultado é que, devido a robustez do método ele se mostra eficiente quando aplicado para o sistema não linear. No método LQR é necessário conhecer todos os estados para a realimentação, caso não seja possı́vel medir todos os estados podemos usar um controlador baseado em observador como é o caso do método LQG. Para o método LQG implentado para o satélite equipado com rodas, considerou-se um modelo com incertezas estocásticas na saı́da e na dinâmica, tornando o estudo mais realista. As medidas dos sensores foram corrompidas com ruı́dos de distribuição gaussiana, a partir dos ajustes do filtro de Kalman obtemos uma boa estimação do estado, como pode se ver nos resultados das simulações. Apesar da perda de robustez que a inserção do filtro pode causar os resultados mostram um desempenho próximo do obtido pelo método LQR, o que evidencia que a seleção das ponderações do controle sugeridas por Kristiansen (2000) é conveniente. Os controladores baseados em energia usando realimentação de atitude e velocidade angular, para o modo de estabilização, mostram-se eficientes, apresentando como vantagem a simplicidade na implementação. A lei de controle utilizando apenas realimentação de velocidade angular não garante, para outras condições iniciais, a estabilização e o controle em três eixos de interesse, sendo necessário informações da atitude. Entretanto ainda é uma alternativa que deve ser considerada, para os casos em que nenhuma informação da atitude é disponı́vel para o sistema de controle. Os resultados mostram que a lei de controle apresentada é uma boa alternativa para o sistema de controle com realimentação de velocidade angular e atitude, obtendo um desempenho razoável para o controle em três eixos. Em ambos os casos o uso de bobinas apresenta uma opção atrativa devido ao baixo consumo de energia. Uma vantagem em relação a metodologia LQR é que mesmo sob manobras de grandes ângulos, a eficiência do sistema de controle de atitude na estabilização do sistema é garantida. Realizou-se o estudo de viabilidade do uso de rodas para manobras de grandes ângulos durante o modo de aquisição. Para esse modo de operação as leis de controle empregadas foram: 1) LQR tracking e 2) controle PD. Para o controle PD foi utilizado a metodologia LQ para o cálculo dos ganhos do controlador. As si143 mulações mostram desempenhos muito próximos dos dois controladores devido ao uso da metodologia LQ para o controlador PD. O método LQR utiliza o estado como realimentação e o PD apenas informações da atitude, entretanto tem a desvantagem de necessitar derivar esses sinais. Os resultados mostram que é viável e factı́vel o uso das rodas especificadas para manobras de aquisição. E mesmo para as manobras de grandes ângulos os controladores que são baseados na planta linear apresentam um bom desempenho. Finalizando esse trabalho contribuiu para a área de dinâmica e controle de atitude no sentido de que: • apresenta um estudo de referenciais e parâmetros mais usados para representar a atitude; • apresenta o modelo matemático da dinâmica de atitude para veı́culos espaciais contendo rodas de reação e bobinas magnéticas; • discute e apresenta as técnicas de controle ótimo (LQR, LQR tracking é LQG) bem como as leis de controle desenvolvindas pelo método de energia; • desenvolve as leis de controle utilizando as técnicas referidas no ı́tem anterior para aplicação de dois conceitos de sistemas de controle para satélites estabilizados em três eixos. Um procedimento que utiliza uma combinação de rodas de reação e bobinas magnéticas como atuadores e outro que utiliza somente bobinas magnéticas com atuadores; • simula em ambiente MATLAB/SIMULINK o controle de satélites estabilizados em três eixos, utilizando diferentes leis de controle para as concepções de sistemas de controle; • disponibiliza um pacote de sofware desenvolvido na plataforma MATLAB/SIMULINK que poderá ser utilizado no futuro para o estudo de dinâmica e controle de atitude. 8.1 Sugestões para Tabalhos Futuros Referente ao modelo seria propõe-se além do modelo do torque de gradiente de gravidade os modelos dos torques de pressão de radiação, arrasto atmosférico, dipolo 144 residual e o efeito jitter, associado a fontes de perturbação interna. Outra sugestão seria a implementação do método LQG para o satélite equipado com bobinas e a apresentação rigorosa das provas de estabilidade para os controladores baseados em energia usando a teoria de Lyapunov. Sugere-se também inserir um filtro para estimação dos dados dos magnetômetros para o projeto do ACS que utiliza bobinas. Para uma simulação mais realista sugere-se: • incluir o modelo da roda de reação; • incluir o modelo do magnetômetro; • testes de robustez através de Monte Carlo. 145 146 REFERÊNCIAS BIBLIOGRÁFICAS Alfriend, T. K. Magnetic attitude control system for dual-spin satellite. AIAA Journal, v. 13, n. 6, p. 817–822. Arantes, J. G.; Fonseca, I. M. da. A comparasion between quaternions and euler angles for satellite atittude dynamics. In: Colóquio Brasileiro de Dinâmica Orbital, 12, 2004. Ubatuba, SP. Anais... Ubatuba, SP: INPE, 2004. Arantes, J. G.; Fonseca, I. M. da. Three-axis attitude dynamics by using torque coils only. In: Colóquio Brasileiro de Dinâmica Orbital, 12, 2004. Ubatuba, SP. Anais... Ubatuba, SP: INPE, 2004. Ardebil. ASRI program. 2004. Disponı́vel <http://www.asri.org.au/ASRI/index.xml>. Acesso em: out 2004. em: Athans, M. The role and use of stochastic linear quadratic gaussian problem in control system design. IEEE Transactions on Automatic Control, v. 6, p. 529–552, 1971. Auer, W. A. A double gimballed momentum wheel for precision three-axis attitude control. In: AGARD Symposium Guidance and Control Panel, 37, 1983. Florence, Italy. Proceedings... Florence, Italy: AGARD, 1983. Bang, M. J. T. H.; Choi, H. D. Large angle attitude control of spacecraft with actuator saturation. Control engineering practice, n. 11, p. 989–997, 2003. Beichman, C.; Neugebauer, G.; Chester, T. IRAS explanatory supplement. 2004. Disponı́vel em: <http://irsa.ipac.caltech.edu/IRASdocs/iras.html>. Acesso em: 10 Aug 2004. Bhat, S. P.; Dham, A. S. Controllability of spacecraft attitude under magnetic actuation. In: Conference on Decision and Control IEEE, 42, 2003. Maui, Hawai, USA. Proceedings... Maui, Hawai, USA: IEEE, 2003. Bishop, R. Modern control system analysis & design using Matlab & Simulink. Mento Park, California: Eddison Wesley, 1997. Brown, C. D. Spacecraft mission design. Reston, Virginia: Education Series AIAA, 1998. 182 p. 147 Brown, R. G. Introduction to randon signals and applied Kalman filtering. New York: John Wiley & Sons, 1997. Bryson, A. E. Control of spacecraft and aircraft. New jersey: Princeton University Press, 1994. 378 p. Buckingham, O. V. A.; Smirnov, G. V. Magnetic torques for momentum desaturation of space station control moment gyros. Journal Spacecraft, v. 9, n. 6, p. 324–330, 1972. Bushenkov, M. Y. O. A. V.; Smirnov, V. G. Attitude stabilization of a satellite by magnetic coils. Acta Astronautica, v. 50, p. 721–728, 2002. Carrara, V. Modelagem das forças e torques atuantes em satélites. 1982. 153 p. (INPE–2454–TDL/094). Dissertação (Mestrado em Ciência Espacial), INPE, São José dos Campos. 1982. Carvalho, H. de C. EQUARS mission analysis. São José dos Campos, INPE, 2003. 40 p. (E2000-TRP-001v00). Chobotov, V. A. Orbital mechanics. Reston, Virginia: Education Series AIAA, 1996. 447 p. Clausen, L. T. The Orsted satellite project. 2004. Disponı́vel em: <http://web.dmi.dk/fsweb/projects/oersted/homepage.html>. Acesso em: jan 2004. Cohen, V. D. Attitude dynamics of an orbiting electromagnet. Journal Spacecraft, v. 11, p. 252–256, 1973. Distefano, A. R. S. J. J.; Willians, I. J. Sistemas de retroação de controle. São Paulo: McGRAW-HILL, 1972. 478 p. Dorato, C. A. P.; Cerone, V. Linear quadratic control an introdution. Englewood Cliffs: New Jersey: Prentice Hall, 1995. 205 p. Doyle, J. C. Guaranteed margins for lqg regulators. IEEE Transactions on Automatic Control, v. 4, n. 23, p. 756–757, 1978. El-Gohary, A. Optimal stabilization of a rigid body motion using rotors system. Applied Mathematics and Computation, v. 136, p. 229–239, 2003. 148 Escobal, P. R. Methods of orbit determination. New York: John Wiley & Sons, 1965. Fauske, K. M. NCUBE attitude control. Norwegian, Trondhein: Departament of Enginneering Cybernetics, NTNU, dec 2002. Fauske, K. M. Attitude stabilization of an underactuated rigid spacecraft. 2003. 55 p. SIV.ING Thesis, (Departament of Engineering Cybernetics). Norwegian University of Technology and Science, Trondheim, Norwegian. Jan 2003. Fichter, M. S. W.; Zentcraf, P. Control design for generalized normal mode operation of bias momentum satellites. Control Engeneering Practice, v. 4, p. 1355–1360, 1996. Flora, A. L. Projeto de um sistema de controle de atitude de um satélite com apêndices flexı́veis pelos métodos LQG/LTR e H-Infinito. 1995. 223 p. Dissertação (Mestrado em Ciência Espacial), INPE, São José dos Campos. Gökçev, P. T. M. C.; Meerkov, S. M. An lqr/lqg theory for system with saturation actuators. Transactions on Automatic Control, v. 46, n. 10, p. 1529–1542, 2001. Grassi, S. V. M.; Moccia, A. Preliminary design of the attitude control system of a microsatellite for earth observation. Space Technological, Biarritz, France, v. 15, p. 223–230, 1995. Gurman, J. B. The SOHO solar cycle mission. 2004. Disponı́vel em: <http://sohowww.nascom.nasa.gov/publications>. Acesso em: sep 2004. Hamzah, N.; Hashida, Y. Tiungsat-1 momentum wheel commissioning. In: World Engineering Congress, 1, 1999. Kuala Lumpur. Proceedings... Kuala Lumpur: World Eng. Cong., 1999. Hanselman, D.; Littlefield, B. Matlab 6 curso completo. São Paulo: Prentice Hall, 2003. 676 p. Hanselman, D. C.; Benjamin, C. K. MATLAB tools for control system analysis and design. Englewood Cliffs: New Jersey: Princeton-Hall, 1994. 149 Hegrenes, O. Attitude control by means of explicit model predictive control, via multi-parametric quadratic programming. 2004. 120 p. Thesis (Engineering Cybernetics), Trondheim Norwegian. 2004. Heidelberg. Computers, displays and space products. 2004. Disponı́vel em: <http://www.teldix.de/>. Acesso em: ago 2004. Hughes, P. C. Spacecraft attitude dyanmics. New York: John Wiley & Sons, 1986. 564 p. Junkins, J. L.; Carrington, C. K. Time optimal magnetic attitude maneuvers. In: IAA/AAS Astrodynamics conference, 11, 1980. USA. Proceedings... USA: IAA/ASS, 1980. Junkins, J. L.; Turner, J. D. Optimal spacecraft rotational maneuvers. New York: Elsevier Science Publishers B. V., 1986. 515 p. Kaplan, M. H. Modern spacecraft dynamic & control. New York: John Wiley & Sons, 1976. 415 p. Khalil, H. K. Nonlinear system. [S.l.: s.n.], 2000. Kim, H. L. B. J.; Choi, S. D. Three axis reaction wheel attitude control system for KITSAT-3 microsatellite. Taejon, Korea: Sattelite Technology Reserch Center, KAIST, 1999. Kirk, D. E. Optimal control theory: an introduction. Englewood Cliffs: New Jersey: Princeton Hall, 1970. Kristiansen. Norwegiann micro salellite. Cybernetics), Trondheim Norwegian. 2000. 2000. Thesis (Engineering Kuga, H. K.; Guedes, U. T. V. Dinâmica de atitude para satélites estabilizados por rotação. São José dos Campos, SP, 1987. (INPE-4403-TVTE/275). Kuga, H. K.; Kondapalli, R. R. Introdução à mecânica orbital. São José dos Campos, 1995. 73 p. (INPE-5615-PUD/064). Kuga L. D., F. H. K.; Guedes, U. T. V. Simulação de atitude de manobras para o satélite brasileiro estabilizado por rotação. São José dos Campos, SP, 1987. (INPE-4271-PRE/1143). 150 Kwakernaak, H.; Sivan, R. Linear optimal control system. New York: John Wiley & Sons, 1972. 564 p. Larson, W. J.; Wertz, J. R. Space mission analysis and design. Torance, California: Space Technology Series, 1992. 865 p. Legg, V. E. Survey of magnetic material and applications in the telephone system. 2003. Disponı́vel em: <http://www.telepsystem.com/reporter/index.xml>. Acesso em: dez 2003. Levine, S. W.; Leonard, N. E. Using Matlab to analyse and control system. New York: Addison Wesley Publishing, 1995. Maciejowski, J. M. Multivariable feedback design. New York: Addison Wesley Publishing, 1989. Macmillan, S.; Quinn, J. M. The derivation of world magnetic model 2000. London: British Geological Survey, 2000. 278 p. (Geomagnetism Series – WM/00/17R). Marteau, P. P. F.; Psiaki, M. Active magnetic control system for gravity gradient stabilized spacecraft. In: Annual AIAA/USU conference for small satellites, 2, 1988. Logan (Utah), USA. Proceedings... Logan (Utah), USA: AIAA, 1988. Marteau, S. B. G. F.; Rogers, E. Attitude determination and control for small spacecraft. In: UKACC International Conference on Control, 13, 1996. U.K. Proceedings... U.K.: IEE, 1996. Martins Neto, A. F. Atitude e seu controle. In: Kuga, A. F. B. de A. P. . H. K. (Ed.). Fundamentos de tecnologia espacial. São José dos Campos, INPE: INPE, 2001. v. 1, p. 65–79. Mclean, S. The world magnetic model and associated software. Washington, D.C., National Oceanic & Atmospheric Administration (NOAA). Disponı́vel em: <http://www.ngdc.noaa.gov>. Acesso em: out 2004. Moore, J. B.; Anderson, D. O. B. Optimal control linear quadratic methods. Englewood Cliffs: New Jersey: Princeton Hall, 1990. 380 p. Moscati, N. R. Projeto de um sistema de controle de atitude (três eixos) de satélites utilizando a metodologia LQG/LTR. 1992. 218 p. 151 (INPE–5473–TDI/504). Dissertação (Mestrado em Ciência Espacial), INPE, São José dos Campos. 1992. Musser, L. K.; Ebert, L. W. Autonomous spacecraft attitude control using magnetic torquing only. In: AIAA Guidance, Navigation, and Control Conference, 12, 1989. NASA Goddard Space Flight Center, Greenbelt. Proceedings... NASA Goddard Space Flight Center, Greenbelt: NASA, 1989. p. 23–38. Ogata, K. Engenharia de controle moderno. Rio de Janeiro: LTC - Livros Técnicos e Cientı́ficos Editora, 1998. 812 p. Overby, E. J. Attitude control for the norwegian student satellite nCube. 2004. 74 p. Thesis (Engineering Cybernetics), Trondheim Norwegian. 2004. Perkel, H. Stabilite - three axis attitude control system utilizing a single reaction wheel. In: AIAA Comunication Satellite System Conference, 1, 1966. Washington, D.C. Proceedings... Washington, D.C.: AIAA, 1966. p. 375–400. Pilchowski, C. C. W. S. H. U.; Ferreira, D. D. L. Introdução à mecânica celeste. São José dos Campos, INPE, 1981. (INPE-COM.4/RPE C.D.U.:521.3). Pilchowski, H. U. Sensores de atuadores. In: Kuga, A. F. B. de A. P. . H. K. (Ed.). Fundamentos de tecnologia espacial. São José dos Campos: INPE, 2001. v. 1, p. 50–64. Psiaki, L. M. Magnetic torquer attitude control via asymptotic periodic linear quadratic regulation. Journal of Guidance, Control, and Dynamics, v. 24, p. 386–394, 2001. Quirelli, I. M. P. Spin stabilized satellite attitude propagation. 2002. Master Thesis, UNESP, Guaratinguetá, SP. 2002. Rodrigues, D. S. S.; Zanardi, M. C. Spacecraft attitude propagation with different representation. In: Kuga, A. F. B. de A. P. . H. K. (Ed.). Advances in space dynamics 4: celestial mechanics and astronautics. São José dos Campos: INPE, 2004. v. 1, p. 143–150. 152 Roma, A. M. Análise dinâmica e controle de um satélite artificial com paı́neis solares flexı́veis. 1991. 177 p. (INPE–5220–TDL/436). Dissertação (Mestrado em Ciência Espacial), São José dos Campos. 1991. Shigehara, M. Geomagnetic attitude control of an axisymemtric spinning satellite. Journal Spacecraft, v. 9, n. 6, p. 623–635, 1972. Silani, E.; Lovera, M. Magnetic spacecraft attitude control survey and some new results. Control Engineering Practice, v. 1, p. 1, 2003. Souza, L. C. G. Controle de atitude de um satélite artificial através da extensão da teoria do regulador linear quadrático. 1987. 60 p. (INPE–4407–TDL/304). Dissertação (Mestrado em Ciência Espacial), INPE, São José dos Campos. 1987. Souza, M. L. O. Estudo e desenvolvimento de um sistema de controle de atitude ativo em três eixos para satélites artificiais usando atuadores pneumáticos a gás frio e volantes a reação. 1981. (INPE–2000–TDL/042). Dissertação (Mestrado em Ciência Espacial), INPE, São José dos Campos. 1981. Souza, P. N. Análise, projeto, construção e testes de um modelo de roda de reação para aplicações espaciais. 1987. 185 p. (INPE–4358–TDL/299). Dissertação (Mestrado em Ciência Espacial), INPE, São José dos Campos. 1987. Spencer, T. M. Automatic magnetic control of a momentum-biased observatory in equatorial orbit. Journal Spacecraft, v. 14, n. 4, p. 211–218, 1977. Spindler, K. Single attitude control laws for momentum-wheel actuated spacecraft. In: International Symposium Space Dynamics, 14, 2000. Biarritz, France. Proceedings... Biarritz, France: AIAA, 2000. Tamura, T. ASCA Measurements of the Gravitational Potential Profile in the Central Region of Galaxy Clusters. Tese (Ph.D Thesis) — Graduate School of Science, University of Tokyo, Tokyo. 1998. Trivelato, G. da C. Controle de rodas de reação através de técnicas digitais usando modelos de referência. 1988. 209 p. (INPE–4618–TDL/335). Dissertação (Mestrado em Ciência Espacial), INPE, São José dos Campos. 1988. 153 Ulrich, H. Satélites Artificiais movimento de atitude. Notas de aula. 2004. Varatharajoo, R.; Fasoulas, S. The combined energy and attitude control system for small satellites earth observation missions. Boston, USA: Sattelite Technology Reserch Center, 1975. Wang, P.; Shtessel, B. Y. B. Satellite attitude control using only magnetictorques. In: AIAA Guidance, Navigation, and Control Conference, 12, 1998. Boston. Proceedings... Boston: AIAA, 1998. Wells, L. S. J. G.; Jeans, T. Canada’s smallest satellite: The canadian advanced nanospace experiment (canx-1). In: Annual AIAA/USU Conference on Small Satellites, 16, 2002. Logan, Utah. Proceedings... Logan, Utah, 2002. Wertz, J. R. Spacecraft attitude determination and control. London, England: D. Reideil Publishing Company, 1978. 861 p. Whitford, C.; Forrest, D. The catsat attitude control system. In: AIAA/USU Conference on Small Satellites, 12, 1998. USA. Proceedings... USA: AIAA, 1998. Wie, B. Space vehicle dynamics and control. Reston, Virginia: AIAA Education Series, 1998. 661 p. Wie, H. W. B.; Arapostathis, A. Quaternion feedback regulator for spacecraft eigenaxis rotaions. Journal of Guidance, Control and Dynamics, v. 12, n. 3, p. 375–380, 1989. Wie, P. M. B. B. Quaternion feedback for spacecraft large angle maneuvers. Journal of Guidance, Control and Dynamics, v. 3, n. 3, p. 360–365, 1985. Wilson, D. Attitude and Orbit Control Using the Spacecraft Control Toolbox v4.6. Princeton - New Jersey, Dec 2000. Wisniewski, R. Satellite attitude control using only electromagnetic actuation. Tese (Ph.D Thesis) — Aalborg University, Departament of Control Engineering, Danish. 1996. Wisniewski, R. Linear time varying approach to satellite control using only electromagnetic actuation. AIAA Guidance, navigation control, v. 11, p. 243–251, 1997. 154 Wisniewski, R.; Blanke, M. Fully magnetic attitude control for spacecraft subject to gravity gradient. Automatica, v. 35, p. 1201–1214, 1999. Wisniewski, R.; Markley, F. L. Optimal magnetic attitude control. In: IAFC World Congress, 14, 1999. USA. Proceedings... USA: IAFC, 1999. p. 313–318. Wright, P. S.; Wong, H. S. An overview of sensors in spacecraft engineering. New York, USA: Sattelite Technology Reserch Center, 1989. Yairi, T. On-board reconfigurable attitude control system with optimazation. In: International Symposium on Space Technology and Science,19,1994. Yokohama, Japan. Proceedings... Yokohama, Japan: ISTS, 1994. Zanardi, M. C.; Assis, S. C. de; Kuga, H. K. Torque residual médio com modelo de quadripolo. In: Congresso Temático de Dinâmica, Controle e Aplicaçòes, 3, 2004. Ilha Solteira - SP. Anais... Ilha Solteira - SP: Série Arquimedes, 2004. Zanardi, M. C.; Quirelli, I. M. P.; Kuga, H. K. Analytical attitude propagation of the spin stabilized earth artificial satellite. In: International Symposium of Space Flight Dynamics, 17, 2003. Moscou - Rússia. Proceedings... Moscou Rússia: CD-ROM, 2003. Zanardi, M. C.; Quirelli, I. M. P.; Kuga, H. K. Torques magnéticos: Aplicações à satélites estabilizados por rotação. In: Congresso Temático de Dinâmica, Controle e Aplicaçòes, 2, 2003. São José dos Campos - SP. Anais... São José dos Campos - SP: Série Arquimedes, 2003. p. 3167–3176. 155 156 APÊNDICE A CÁLCULO DO GRADIENTE DE GRAVIDADE No desenvolvimento do modelo do torque devido ao gradiente de gravidade consideram-se as seguintes aproximações, válidas para satélites em órbita da Terra em sua maioria. • somente o campo gravitacional da Terra é considerado; • a Terra é considerada esférica com distribuição de massa uniforme; • as dimensões do satélite são pequenas em comparação com a distância Terra - satélite; • o satélite é um corpo simples . A partir das aproximações 1, 2 e 3 o torque devido ao gradiente de gravidade pode ser expresso por Z τ g = −µ sat r×R dm R3 (A.1) onde µ é a constante geo-gravitacional, R é o vetor distância do centro da Terra ao elemento de massa dm do satélite e r é o vetor distância do centro de massa do satélite ao elemento de massa dm. Tem-se que R = |Rs + r| (A.2) onde Rs é o vetor distância do centro de massa do satélite ao centro de massa da Terra. Substituindo a Equação (A.2) em (A.1), expandindo o termo |Rs + r|−3 em série 2 |r| de Taylor e desprezando os termos de ordem igual e superior a dois, O |R , s| |r| admissı́vel com base na aproximação 4 ( |R << 1) obtemos s| 157 µ τg = 5 Rs Z (Rs · r) (r × Rs ) dm (A.3) sat Manipulando a Equação (A.3) como segue (Wie, 1998) τg Z µ = −3 5 Rs × r (r · Rs ) dm Rs Zsat µ = −3 5 Rs × rrdm · Rs Rs sat (A.4) (A.5) Usando a notação dyadic, pormenorizadas em Hughes (1986) e Wie (1998), o tensor de inércia pode ser escrito como Z J= r2 Î − rr (A.6) sat onde Î = x̂x̂ + ŷ ŷ + ẑ ẑ é o dyadic unitário, e x̂, ŷ e ẑ são os vetores unitários/versores que formam a base do sistema do satélite. Substituindo (A.6) na Equação (A.5), temos τg Z µ 2 = −3 5 Rs × r Î − J · Rs Rs sat Z µ µ = −3 5 Rs × r2 Îdm · Rs + 3 5 Rs × J · Rs Rs Rs sat (A.7) (A.8) de (A.6) e da relação Î · Rs = Rs obtemos τg = 3 µ Rs × J · Rs Rs5 158 (A.9) O torque devido ao gradiente de gravidade pode ser expresso por Rs µ Rs × J · Rs3 Rs Rs µ = 3 3 R̂s × J · R̂s Rs τg = 3 (A.10) (A.11) onde R̂s é o vetor unitário na direção Rs . Escrevendo o torque de gradiente de gravidade no referencial orbital ou VLHL, resulta τ og = 3 µ ẑo × J · ẑo Rs3 (A.12) onde ẑo é o vetor intário na direção nadir/vertical local. A velocidade orbital média é dada por r ωo = µ Rs3 (A.13) Substituindo na Equação (A.12) obtemos τ og = 3ωo2 ẑo × J · ẑo (A.14) Usando a notação do operador anti-simétrico para o produto vetorial, reescrevemos a Equação (A.14) τ og = 3ωo2 S (ẑo ) J · ẑo Expressando o vetor unitário ẑo no referencial do satélite (x̂, ŷ, ẑ) temos 159 (A.15) ẑo = cb3x x̂ + cb3y ŷ + cb3z ẑ (A.16) onde cb3i , (i = x, y, z) são os cosenos diretores do vetor initário ẑo nas direções x̂, ŷ, ẑ. Desenvolvendo a Equação (A.15) b c3x 0 −cb3z cb3y b b 2 b b τ g = 3ωo c3z 0 −c3x J c3y −cb3y cb3x 0 cb3z (A.17) Como os eixos do sistema de referência do satélite coincide com os eixo principais de inércia, temos que os produtos de inércia são zeros, temos que a Equação (A.17) resulta b c3x Jx 0 −cb3z cb3y b b 2 b b τ g = 3ωo c3z 0 −c3x c3y Jy b b 0 cb3z Jz −c3y c3x (A.18) Desenvolvendo a Equação (A.18) obtemos o modelo do torque de gradiente de gravidade expresso no referencial do satélite (Jz − Jy ) cb3y cb3z τ bg = 3ωo2 (Jx − Jz ) cb3x cb3z (Jy − Jx ) cb3x cb3y (A.19) Podemos escrever o gradiente de gravidade em termos dos ângulos de Euler/de atitude; roll (ψ), pitch (θ) e yaw (φ) como segue Para a seqüência de rotação 3 − 2 − 1 os cosenos diretores cb3 são dados por (Wertz, 1978) 160 cb3x = − sin θ (A.20) cb3y = sin φ cos ψ (A.21) cb3z = cos φ cos ψ (A.22) Substituindo na Equação (A.19) obtemos o modelo do torque de gradiente de gravidade em função dos ângulos de Euler (Jz − Jy ) sin φ cos φ cos θ2 τ bg = 3ωo2 (Jx − Jz ) cos φ cos θ sin θ (Jy − Jx ) sin φ sin θ cos θ (A.23) Usando a Equação (4.11) do Capı́tulo (4), Seção (4.2.2), podemos expressar o gradiente de gravidade em função dos parâmetos de Euler/quaternions como (Jz − Jy ) (1 η + 2 3 ) (1 − 221 − 222 ) τ bg = 6ωo2 (Jx − Jz ) (1 3 − 2 η) (1 − 221 − 222 ) 2 (Jy − Jx ) (1 3 − 2 η) (1 η + 2 3 ) (A.24) Linearização da Expressão do Torque de Gradiente de Gravidade No modo de operação onde são realizadas manobras de pequenos ângulos, é razoável aproximar as equações por equações lineares, válidas em torno de um ponto de operação. Considerando θ, φ e ψ menores que até 15 graus, podemos fazer as seguintes aproximações sin φ ≡ φ, sin θ ≡ θ, cos θ ≡ 1 e θφ ≡ 0. Fazendo essas aproximações, o modelo do torque de gradiente de gravidade, dado pela Equação (A.23) é aproximado por (Wie, 1998, Kaplan, 1976) (Jz − Jy ) φ τ bg = 3ωo2 (Jx − Jz ) θ 0 161 (A.25) Em termos dos parâmetros de Euler o modelo do torque de gradiente de gravidade pode ser aproximado por (Jz − Jy ) 21 τ bg = 3ωo2 (Jx − Jz ) 22 0 162 (A.26) APÊNDICE B PROPAGAÇÃO DA ÓRBITA E TRANSFORMAÇÕES Esse apêndice descreve de maneira detalhada a teoria necessária para o cálculo da propagação da órbita Kepleriana (modelo de dois corpos). São descritas ainda as transformações entre os diferentes referênciais envolvidos e o cálculo da longitude, latitude e altura do ponto sub-sátelite. Todas as transformações e o cálculo da órbita são implementados em MATLAB. Para a propagação do campo magnético terrestre (modelo IGRF) é necessário incluir: • cálculo/propagação da órbita (Kepleriana); • matriz de transformação do sistema inercial para o sistema do satélite; • obtenção da latitude, longitude e a altura do ponto sub-satélite; A partir desses cálculos podemos obter as componentes do campo magnético terrestre, espressas no referencial do satélite. Nota: O modelo IGRF, disponı́vel nas linguagens C e FORTRAN foi convertido para ambiente MATLAB e SIMULINK, usando S-functions. Detalhes desse procedimento é encotrado na documentação da mathworks, disponı́vel em www.mathworks.com. B.1 B.1.1 Cálculo da Órbita Posicionamento de Satélites - Problema Direto O problemo direto consiste em dados os elementos Keplerianos (a, e, i, ω, Ω, M ) determinar à posição do satélite em coordenadas cartesianas X = [Xi Zi Zi ]T (Kuga e Kondapalli, 1995). Onde X é o vetor posição no referencial inercial (ECI). 163 Figura B.1- Referenciais (inercial, da órbita) e os elementos Keplerianos (i, ω, Ω). Pela Figura (B.1) pode-se notar que os elementos Keplerianos que representam os ângulos de Euler são: • Ω é a ascensão do nodo ascendente 0o ≤ Ω ≤ 360o ; • i é a inclinação da órbita em relação ao equador −90o ≤ i ≤ 90o ; • ω é o argumento do perigeu 0o ≤ ω ≤ 360o . Esses elementos Keplerianos definem a orientação da órbita. Os elementos Keplerianos (a, e) definem o tamanho e tipo (elı́ptica, circular, hiperbólica) da órbita. Sistema da órbita (OOF): O sistema da órbita definido (Xo , Yo , Zo ) é um sistema de coordenadas com origem no centro de massa da Terra. O eixo Xo aponta na direção do perigeu Π, o eixo Zo aponta na direção normal ao plano da órbita. O eixo Yo é obtido pela regra da mão direita. B.1.2 Equação de Kepler As coordenadas do satélite em relação ao sistema da órbita são dadas por 164 Xo = a (cos u − e) √ Yo = a 1 − e2 (sin u) (B.1) Zo = 0 (B.3) (B.2) onde a é o semi-eixo maior, u é a anomalia excêntrica e e é a excentricidade da órbita. A Equação de Kepler é dada por (Brown, 1998) M = u − e sin u (B.4) onde M é a anomalia média dada por M = n (t − to ) (B.5) n é o movimento médio, que é dado por r µ a3 (B.6) µ = GM (B.7) n= e onde G é a constante gravitacional e M a massa do corpo central. Para a Terra 3 temos µ = 398600.4 km s2 165 B.1.3 Matriz de Rotação Dados os elementos orbitais (i, Ω, ω) temos que a matriz de transformação/rotação que fornece as coordenadas do sistema inercial em função das coordenadas do sistema da órbita é dada por (Kuga e Kondapalli, 1995). R (i, Ω, ω) = RZo (−Ω) RXo (−i) RZo (−ω) (B.8) Sendo a relação/transformação dos sistemas expressa por Xi Xo Yi = R (i, Ω, ω) Yo Zi Zo (B.9) Temos, para cada seqüência de rotações, as matrizes cos ω − sin ω 0 RZo (−ω) = sin ω cos ω 0 0 0 1 (B.10) cos Ω − sin Ω 0 RZo (−Ω) = sin Ω cos Ω 0 0 0 1 (B.11) 1 0 0 RXo (−i) = 0 cos i − sin i 0 sin i cos i (B.12) Substituindo na Equação (B.8) obtemos. 166 cos Ω cos ω − sin Ω cos i sin ω − cos Ω sin ω − sin Ω cos i cos ω sin Ω sin i R (i, Ω, ω) = sin Ω cos ω + cos Ω cos i sin ω − sin Ω sin ω + cos Ω cos i cos ω − cos Ω sin i sin i sin ω sin i cos ω cos i (B.13) Portanto, de (B.1), (B.2), (B.3) e (B.13) em (B.9), obtemos as coordendas do satélite no sistema inercial. Xi (cos Ω cos ω − sin Ω cos i sin ω) Xo − (cos Ω sin ω − sin Ω cos i cos ω) Yo Yi = (sin Ω cos ω + cos Ω cos i sin ω) Xo − (sin Ω sin ω + cos Ω cos i cos ω) Yo Zi sin i sin ωXo + sin i cos ωYo (B.14) Órbita Circular Para uma órbita circular temos que: a = R⊗ + h e ∼ = 0 (B.15) ω = 0 (B.17) (B.16) Note que a referência para o inı́cio do cálculo da órbita circular é tomado sobre o nodo. R⊗ é o raio médio da Terra e h a altitude do satélite. Da Equação de Kepler (B.4) resulta. M =u=f logo de (B.5) temos 167 (B.18) u = n (t − to ) (B.19) onde n = ωo = velocidade orbital constante. As equações (B.1), (B.2) e (B.3) se reduzem a Xo = a cos u (B.20) Yo = a sin u (B.21) Zo = 0 (B.22) e a matriz de transformação (B.13) se reduz a cos Ω − sin Ω cos i sin Ω sin i R (i, Ω, ω) = sin Ω cos Ω cos i − cos Ω sin i 0 sin i cos i (B.23) Substituindo na Equação (B.9), obtemos as coordenadas dos satélite no sistema inercial para o caso de uma órbita circular. Xi a cos Ω cos u − a sin Ω cos i sin u Yi = a sin Ω cos u + a cos Ω cos i sin u Zi a sin i sin u (B.24) A Tabela (B.1) apresenta o programa em MATLAB para a propagação da órbita circular. 168 Tabela B.1- Programa MATLAB para o cálculo da órbita circular. function [X] = SimOrbitaECI(par) % by Gilberto Arantes 2004 INPE % % [X] = SimOrbitaECI(par) Propaga uma orbita circular % % Calculo da posicao do satelite no referencial ECI % % Input % par = [ h i asc nodo u] % h = altura do satelite (km) % i = inclinacao da orbita (graus) % asc nodo = ascensao do nodo ascendente (graus) % u = anomalia excentrica (rad) % % Output % X = vetor posicao no referencial inercial (km) % h = par(1); i = par(2); asc nodo = par(3); u = par(4); r = 6378 + h; % a = distancia do centro da Terra ao satelite em km % Calculo da matriz de tranformacao do sistema orbital para o inercial (ref.: Kuga, 1995) fator = pi/180; a = 0; % Orbita circular: argumento do perigeu e zero b = fator*i; c = fator*asc nodo; R = [(cos(c)*cos(a)-sin(c)*cos(b)*sin(a)) (-cos(c)*sin(a)-sin(c)*cos(b)*cos(a)) (sin(c)*sin(b)) (sin(c)*cos(a)+cos(c)*cos(b)*sin(a)) (-sin(c)*sin(a)+cos(c)*cos(b)*cos(a)) (-cos(c)*sin(b)) sin(b)*sin(a) sin(b)*cos(a) cos(b) ]; % Posicao no referencial orbital xo = r*cos(u); yo = r*sin(u); zo = 0; Xo=[xo;yo;zo]; % Coordenadas no Sistema da Orbita % Posicao no referencial inercial ECI X = R*Xo; % X = X [X Y Z]’ % end B.2 Matriz de Rotação do Sistema Inercial - Sistema do Sátelite O cálculo da matriz de transformação do sistema inercial ECI (definido em 4.1.1) para o sistema do satélite é mostrado nessa Seção. As etapas do procedimento para o cálculo da matriz de transformação são: • obtemos a matriz de transformação do referencial inercial (ECI) para o pseudo 169 OF orbital (POF) (RPECI ); • obtemos a matriz de transformação do referencial pseudo orbital (POF) para o orbital (OF)(ROF P OF ); • obtemos a matriz de transformação do referencial orbital (OF) para o sistema do satélite (BF) (RBF OF ). Obtida cada uma das matrizes de rotação, temos que a matriz de transformação do sistema inercial para o sistema do satélite será dada por BF OF P OF RBF ECI = ROF RP OF RECI B.2.1 (B.25) OF Matriz RPECI Referencial pseudo orbital: O referencial pseudo orbital (POF) definido por 0 0 0 Xo , Yo , Zo é um sistema de coordenadas com origem no centro da Terra. O eixo 0 0 Xo aponta na direção do centro de massa do satélite, o eixo Zo aponta na direção 0 normal a plano da órbita. O eixo Yo é obtido pela regra da mão direita (Ulrich, 2004). Note que o referencial pseudo orbital é um referencial girante, assim como o referencial orbital, definido em 4.1.2. A Figura (B.2) ilustra o referido referencial. Figura B.2- Referencial inercial (ECI) e pseudo orbital (POF). 170 A partir da Figura (B.2) podemos ver que a seqüência de rotações que leva o sistema inercial para o sistema pseudo orbital é dada por RZi (u) ← RXi (i) ← RZi (Ω), ou seja, primeira rotação de um ângulo Ω em torno do eixo Zi , segunda rotação de um ângulo i em torno do eixo Xi e terceira rotação de um ângulo u em torno do eixo Zi , essa seqüencia é válida para uma órbita circular. Temos para cada seqüência de rotação cos u sin u 0 RZi (u) = − sin u cos u 0 0 0 1 (B.26) cos Ω sin Ω 0 RZi (Ω) = − sin Ω cos Ω 0 0 0 1 (B.27) 1 0 0 RXi (i) = 0 cos i sin i 0 − sin i cos i (B.28) Temos que a matriz de transformação é dada por R (i, Ω, u) = RZi (u) RXi (i) RZi (Ω) (B.29) 0 Xo Xi 0 Yo = R (i, Ω, u) Yi 0 Zo Zi (B.30) onde Substituindo a Equação (B.26), (B.27) e (B.28) em (B.29), obtemos a matriz de 171 rotação que leva o sistema inercial (ECI) para o sistema pseudo orbital (POF) OF RPECI B.2.2 cos u cos Ω − cos i sin u sin Ω cos Ω cos i sin u + cos u cos Ω sin i sin u = − cos Ω sin u − cos i cos u cos Ω cos Ω cos i cos u − sin u sin Ω sin i cos u sin i sin Ω − cos Ω sin i cos i (B.31) Matriz ROF P OF O referencial orbital foi definido em (4.1.2), a Figura (B.3) ilustra a orientação do referencial orbital em relação ao referencial pseudo orbital, note que o referencial pseudo orbital é representado no centro de massa do satélite. Figura B.3- Orientação do referencial orbital (OF) e do referencial pseudo orbital (POF). A seqüência de rotações que descreve o referencial pseudo orbital em relação ao referencial orbital é RXo0 (−90o ) ← RZo0 (90o ), onde 0 1 0 RZo0 (90o ) = −1 0 0 0 0 1 172 (B.32) 1 0 0 RXo0 (−90o ) = 0 0 −1 0 1 0 (B.33) Temos que a matriz de transformação é dada por o o ROF P OF = RXo0 (−90 ) RZo0 (90 ) (B.34) onde xo OF yo = RP OF zo 0 Xo 0 Yo 0 Zo (B.35) Substituindo as equações (B.32) e (B.33) em (B.34), obtemos a matriz de transformação do referencial pseudo orbital para o referencial orbital ROF P OF B.2.3 0 1 0 = 0 0 −1 −1 0 0 (B.36) Matriz RBF OF A matriz de rotação do sistema orbital para o sistema do corpo é obtida a partir dos ângulos de atitude ou ângulos de Euler. A seqüência de rotações escolhida foi 3 − 2 − 1, cuja matriz de atitude é dada por (Wertz, 1978) 173 RBF OF cos θ cos ψ cos θ sin ψ − sin θ = − cos φ sin ψ + sin φ sin θ cos ψ cos ψ cos φ + sin ψ sin θ sin φ sin φ cos θ sin ψ sin φ + cos φ sin θ cos ψ − sin φ cos ψ + cos φ sin θ sin ψ cos φ cos θ (B.37) Matriz RBF ECI Substituindo as equações (B.31), (B.36) e (B.37) na Equação (B.25), e com o auxı́lio do manipulador simbólico do programa MATHEMATICA, obtemos finalmente a matriz de transformação do sitema inercial para o sistema do satélite. RBF ECI c11 c12 c13 = c21 c22 c23 c31 c32 c33 (B.38) com c11 = −cθsisψsΩ − cθcψ (cΩsu − cicusΩ) − sθ (cisusΩ − cucΩ) c12 = cθsisψcΩ + cθcψ (−sΩsu − cicucΩ) − sθ (−cisucΩ − cusΩ) c13 = cucθsψsi + sisusθ − cicθsψ (B.39) que são os cosenos diretores de Xi em relação aos eixos do satélite (x, y, z) c21 = −sisΩ (cφ + sθsφsψ) + (cψsθsφ − cφsψ) (−cΩsu − cicusΩ) + cθsφ (−cucΩ + cisusΩ) c22 = sicΩ (cφ + sθsφsψ) + (cψsθsφ − cφsψ) (−sΩsu + cicucΩ) + cθsφ (−cusΩ − cisucΩ) c23 = −cθsisusφ + cusi (cψsθsφ − cφsψ) − ci (cφcψ + sθsφsψ) (B.40) que são os cosenos diretores de Yi em relação aos eixos do satélite (x, y, z) 174 c31 = −sisΩ (−cψsφ + cφsθsψ) + (cφcψsθ + sφsψ) (−cΩsu − cicusΩ) + cθcφ (cisusΩ − cucΩ) c22 = sicΩ (−cψsφ + cφsθsψ) + (cφcψsθ + sφsψ) (−sΩsu + cicucΩ) + cθcφ (−cisucΩ − cusΩ) c23 = −cθcφsisu − ci (cφsθsψ − cψsφ) + cusi (cφcψsθ + sφsψ) (B.41) que são os cosenos diretores de Zi em relação aos eixos do satélite (x, y, z) A matriz RBF ECI é usada para expressar o vetor campo magnético, fornecido em coordenadas inerciais, no sistema do satélite. As Tabelas (B.2),(B.3) e (B.4) apresentam o programa em MATLAB para o cálculo BF da matriz RECI . B.3 Cálculo da Latitude, Longitude e Altura Para o cálculo da latitude, longitude e altura, primeiro expressamos a posição do satélite no referencial cartesiano geocêntrico. Referencial cartesiano terrestre ou geocêntrico: (FG) (Xg , Yg , Zg ) é um sistema de coordenadas com origem no centro de massa da Terra. O eixo Xo aponta na direção do meridiano de Greenwich, o eixo Zo aponta na direção normal ao plano do equador. O eixo Yo é obtido pela regra da mão direita. A Figura (B.4) ilustra o referencial inercial e o referencial cartesiano geocêntrico (Kuga e Kondapalli, 1995). Na Figura (B.4) θg é o tempo sideral de Greeenwich. O cálculo do tempo sideral de Greenwich é mostrado em Escobal (1965). A matriz de transformação do sistema inercial para o sistema cartesiano geocêntrico é dada por (Pilchowski e Ferreira, 1981) cos θ sin θ 0 G RFECI (θg ) = − sin θ cos θ 0 0 0 1 175 (B.42) BF . Tabela B.2- Programa MATLAB para o cálculo de c1i , (i = 1, 2, 3) da matriz RECI function [R1] = ECI2BF X(par) % by Gilberto Arantes 2004 INPE % % [R1] = ECI2BF X(par) Cosenos diretores X.x X.y X.z % % Calculo dos cosenos diretores % % Input % par = [orb atitude] % % atitude [phi teta psi] = angulos de atitude (rad) % orb = [u i asc nodo] = elementos orbitais/Keplerianos (rad, graus, graus] % % Output % R1 = [r11,r21,r31] % r11 = cos(X.x) % r21 = cos(X.y) % r31 = cos(X.z) orb(1) = par(1); %u orb(2) = par(2); %i orb(3) = par(3); %asc nodo % f = pi/180; %graus->radianos at(1) = par(4); %rad phi at(2) = par(5); %rad teta at(3) = par(6); %rad psi % % Calculo dos cosenos diretores % % cos(X.x) r11 = -cos(at(2))*sin(orb(2)*f)*sin(at(3))*sin(orb(3)*f)+... cos(at(2))*cos(at(3))*(-cos(orb(3)*f)*sin(orb(1))-cos(orb(2)*f)*cos(orb(1))*sin(orb(3)*f))-... sin(at(2))*(-cos(orb(1))*cos(orb(3)*f)+cos(orb(2)*f)*sin(orb(1))*sin(orb(3)*f)); % cos(X.y) r21 = cos(at(2))*sin(orb(2)*f)*sin(at(3))*cos(orb(3)*f)+... cos(at(2))*cos(at(3))*(-sin(orb(3)*f)*sin(orb(1))+cos(orb(2)*f)*cos(orb(1))*cos(orb(3)*f)-... sin(at(2))*(-cos(orb(1))*sin(orb(3)*f)-cos(orb(2)*f)*sin(orb(1))*cos(orb(3)*f)); % cos(X.z) r31 = cos(at(2))*sin(orb(2)*f)*cos(at(3))*cos(orb(1))+... sin(orb(2)*f)*sin(orb(1))*sin(at(2))-cos(orb(2)*f)*cos(at(2))*sin(at(3)); % R1 = [r11 r21 r31]; % % end Obtido as coordenadas do satélite no referencial cartesiano geocêntrico podemos obter a latitude φ a longitude λ e a altura h do satélite, como segue. 176 BF . Tabela B.3- Programa MATLAB para o cálculo de c2i , (i = 1, 2, 3) da matriz RECI function [R2] = ECI2BF Y(par) % by Gilberto Arantes 2004 INPE % % [R2] = ECI2BF Y(par) Cosenos diretores Y.x Y.y Y.z % % Calculo dos cosenos diretores % % Input % par = [orb atitude] % atitude [phi teta psi] = angulos de atitude (rad) % obr = [u i asc nodo] = elementos orbitais/Keplerianos (rad, graus, graus] % % Output % R1 = [r12,r22,r32] % r12 = cos(Y.x) % r22 = cos(Y.y) % r32 = cos(Y.z) % f = pi/180; %graus->radianos u = par(1); % rad i = par(2); % graus nodo = par(3); % graus % phi = par(4); %rad phi teta = par(5); %rad teta psi = par(6); %rad psi % % Calculo dos cosenos diretores % % cos(Y.x) r12 = -sin(i*f)*sin(nodo*f)*(cos(phi)*cos(psi)+sin(teta)*sin(phi)*sin(psi))... +(cos(psi)*sin(teta)*sin(phi)-cos(phi)*sin(psi))*(-cos(nodo*f)*sin(u)-cos(i*f)*cos(u)*sin(nodo*f))+... cos(teta)*sin(phi)*(-cos(u)*cos(nodo*f)+cos(i*f)*sin(u)*sin(nodo*f)); % cos(Y.y) r22 = cos(nodo*f)*sin(i*f)*(cos(phi)*cos(psi)+sin(teta)*sin(phi)*sin(psi))+... cos(teta)*sin(phi)*(-cos(u)*sin(nodo*f)-cos(i*f)*sin(u)*cos(nodo*f))+... (cos(psi)*sin(teta)*sin(phi)-cos(phi)*sin(psi))*(-sin(nodo*f)*sin(i*f)+... cos(i*f)*cos(u)*cos(nodo*f)); % cos(Y.z) r32 = -cos(teta)*sin(i*f)*sin(u)*sin(phi)+... cos(u)*sin(i*f)*(cos(psi)*sin(teta)*sin(phi)-cos(phi)*sin(psi))-... cos(i*f)*(cos(phi)*cos(psi)+sin(teta)*sin(phi)*sin(psi)); % R2 = [r12 r22 r32]; % % end Observando a Figura (B.5), temos que 177 BF . Tabela B.4- Programa MATLAB para o cálculo de c3i , (i = 1, 2, 3) da matriz RECI function [R3] = ECI2BF Z(par) % by Gilberto Arantes 2004 INPE % % [R3] = ECI2BF Z(par) Cosenos diretores Z.x Z.y Z.z % % Calculo dos cosenos diretores % % Input % par = [orb atitude] % atitude [phi teta psi] = angulos de atitude (rad) % obr = [u i asc nodo] = elementos orbitais/Keplerianos (rad, graus, graus] % % Output % R3 = [r13,r23,r33] % r13 = cos(Z.x) % r23 = cos(Z.y) % r33 = cos(Z.z) % f = pi/180; %graus->radianos u = par(1); % rad i = par(2); %rad nodo = par(3); %rad % phi = par(4); %rad phi teta = par(5); %rad teta psi = par(6); %rad psi % % Calculo dos cosenos diretores % % cos(Z.x) r13 = -sin(i*f)*sin(nodo*f)*(-sin(phi)*cos(psi)+sin(teta)*cos(phi)*sin(psi))+... (cos(psi)*sin(teta)*cos(phi)+sin(phi)*sin(psi))*(-cos(nodo*f)*sin(u)-cos(i*f)*cos(u)*sin(nodo*f))+... cos(teta)*cos(phi)*(-cos(u)*cos(nodo*f)+cos(i*f)*sin(u)*sin(nodo*f)); % cos(Z.y) r23 = cos(nodo*f)*sin(i*f)*(-sin(phi)*cos(psi)+sin(teta)*cos(phi)*sin(psi))+... cos(teta)*cos(phi)*(-cos(u)*sin(nodo*f)-cos(i*f)*sin(u)*cos(nodo*f))+... (cos(psi)*sin(teta)*cos(phi)+sin(phi)*sin(psi))*(-sin(nodo*f)*sin(u)+cos(i*f)*cos(u)*cos(nodo*f)); % cos(Z.z) r33 = -cos(teta)*sin(i*f)*sin(u)*cos(phi)+... cos(u)*sin(i*f)*(cos(psi)*sin(teta)*cos(phi)+sin(phi)*sin(psi))-... cos(i*f)*(-sin(phi)*cos(psi)+sin(teta)*cos(phi)*sin(psi)); % R3 = [r13 r23 r33]; % % end Xg cos φ cos λ Yg = (R⊗ + h) cos φ sin λ Zg sin φ 178 (B.43) Figura B.4- Referencial inercial e cartesiano geocêntrico. Figura B.5- Longitude (λ), latitude (φ) e altura (h) do satélite. Logo, a longitude será dada por: λ = arctan Yg Xg onde 0 ≤ λ ≤ 2π (B.44) A latitude será dada por: φ = arcsin Zg R⊗ + h onde 179 − π π ≤φ≤ 2 2 (B.45) e a altura dada por: h= q Xg2 + Yg2 + Zg2 − R⊗ (B.46) A Tabela (B.5) apresenta o programa em MATLAB para o cálculo da longitude, latitude do ponto sub-satélite e altura do satélite. 180 Tabela B.5- Programa MATLAB para o cálculo da longitude, latitude do ponto sub-satélite e altura do satélite. function [LLA] = ECI2LLA(par) % by Gilberto Arantes 2004 INPE % % [LLA] = ECI2LAA(par) para uma orbita circular % % Calculo da latitude, longitude do ponto sub-satelite e altitude do satelite % % Input % par = [ teta g X ] % teta g = Tempo sideral em graus % X = [X Y Z] coordenadas no referencial inercial em km % % Output % LLA = [latitude longitude altitude] % latitude em rad % Longitude em rad % altitude em km % teta g = par(1); X(1) = par(2)*1000; % transformando para metros X(2) = par(3)*1000; X(3) = par(4)*1000; % Raio t = 6378000; % R = Raio da Terra em m % % Calculo da matriz de tranformacao do sistema inercial para o terresrte (ref.: Silva, 1981) % fator = pi/180; a = teta g*fator; % R = [ cos(a) sin(a) 0 -sin(a) cos(a) 0 0 0 1]; % Posicao no referencial terrestre x g = R*transpose(X); % x g = [x g y g z g]’ % Calculo da altura h = ((sqrt( x g(1)^ 2 + x g(2)^ 2 + x g(3)^ 2 ) - Raio t))/1000; % Calculo da latitude lat = asin(x g(3)/(Raio t+h)); % Sul (-) Norte (+) % Calculo da longitude long = atan2(x g(2),x g(1)); % Leste (+) Oeste (-) % conversao da saida (0 <= long <= 360) if sign(long) < 0; long = 2*pi+long; else sign(long) >= 0; long = long; end % LLA = [lat long h]; % % end 181 182 APÊNDICE C IMPLEMENTAÇÃO EM SIMULINK Os modelos apresentados no Capı́tulo (5) e as leis de controle apresentadas no Capı́tulo (6) foram implementedas em SIMULINK. O sistema foi dividido nos blocos de: 1) dinâmica que contêm o modo completo do satélite equipado com bobinas e rodas; 2) ambiente que contem o modelo do gradiente de gravidade e do campo magnético; 3) órbita com o cálculo da órbita (modelo kepleriano) do satélite e 3) com os controladores avaliados para os três modos de operação. A implentação é feita para os três modos de operação: • Detumble; • Estabilização; • aquisição. A Figura (C.1) mostra a implemtação em SIMLINK do modo de detumble. A Figura (C.2) mostra o controlador usado para o detumbling. Figura C.1- Implemetação em SIMULINK do Modo de detumble. 183 Figura C.2- Lei de controle Bdot. A Figura (C.3) mostra o modo de estabilização, para o procedimento de controle caracterizado por bobinas, os diferentes controladores, LQR e baseados em energia, que são implentados usando-se multiport switches, Figura (C.4), o que torna fácil selecionar a lei de controle para o modo de estabilização. Figura C.3- Implementação em SIMULINK do modo de estabilização: satélite equipado com bobinas. As Figuras (C.5) e (C.6) mostram o modo de estabilização, com comtrole via rodas de reação. A primeira utilizando o método LQR e a segunda utilizando o método LQG. São mostrados o modelo linear e o completo (não linear) da dinâmica, para efeito de comparação dos resultados relativos à aplicação das leis de controle para ambos os modelos. 184 Figura C.4- Controle LQR e controladores baseados em energia. Figura C.5- Modo de estabilização: satélite equipado com rodas, usando o método LQR. 185 Figura C.6- Modo de estabilização: satélite equipado com rodas, usando o método LQG. A Figura (C.7) mostra o sistema de controle para o modo de aquisição utilizandose um controlador PD. A Figura (C.8) mostra a lei de controle do método LQR tracking implementada para a fase de aquisição. Figura C.7- Modo de aquisição usando o controlador PD. 186 Figura C.8- Lei de controle LQR tracking. Figura C.9- Lei de controle PD. As Figuras (C.10) e (C.11) mostram a implementação em SIMULINK dos modelos: satélite equipado com bobinas e satélite equipado com rodas, respectivamente. Figura C.10- Modelo dinâmico do satélite equipado com bobinas. A Figura (C.12) mostra o modelo das rodas de reação. Note que as saturações são incluı́das. A Figura (C.13) mostra o modelo das bobinas magnéticas. 187 Figura C.11- Modelo dinâmico do satélite equipado com rodas. Figura C.12- Modelo das rodas. A Figura (C.14) mostra o modelo do gradiente de gravidade e a Figura (C.15) mostra o modelo do campo magnético. As Figuras (C.16), (C.17) e C.18) mostram a implemetação da propagação da órbita, do cálculo da latitude e longitude do ponto sub-satélite e das transformações dos referenciais envolvidos, respectivamente. Note que a implemtentação para órbita e as transformações utilizam os arquivos mostrados em MATLAB mostrados na Seção (B). 188 Figura C.13- Modelo das bobinas. Figura C.14- Modelo do torque de gradiente de gravidade. 189 Figura C.15- Modelo do campo magnético. Figura C.16- Propagação da órbita. 190 Figura C.17- Cálculo da latitude, longitude e altura. Figura C.18- Matriz de transformação entre os sistemas de referência ECI e BF. 191 192 APÊNDICE D PROGRAMAS EM MATLAB A seguir são apresentados os arquivos de inicialização com o projeto LQR para os modelos do satélite equipado com rodas e bobinas e o projeto LQG. D.1 Projeto LQR As Tabelas (D.1) e (D.3) mostram os programas em MATLAB para o projeto LQR dos dois modelos analisados. D.2 Projeto LQG A Tabela (D.4) mostra o programa em MATLAB para o projeto LQG, usado no modelo do satélite equipado com rodas. 193 Tabela D.1- Programa MATLAB para o projeto LQR do modelo satélite equipado com rodas. % Por Gilberto Arantes @INPE 2003 % % Gera as matrizes para o projeto do sistema de controle LQR e LQR rastreio % Calculo da equacao matricial de Ricatti e implementacao do LQR % Matrizes A, B, C e D do sistema: satélite equipado com rodas % d/dt(x) = Ax + Bu % y = Cx +du % Inicializacao clc clear all close all % deg2rad = pi/180; rad2deg = 180/pi; mi=3.986e14; % moviemnto medio em (m) r=(6378+750)1000; % raio da orbita circular em (m) n=sqrt(mir^ 3)); % movimento orbital medio em (rad/s) H o = n*0.01; % Momento angular nominal kgm2/s % Momentos e produtos de inercia (14/10/2003 de Sebastiao) % Ixx = 13.31 kg.m2 Ixy = 0.37 kg.m2 % Iyy = 14.22 kg.m2 Ixz = 0.39 kg.m2 % Izz = 11.2 kg.m2 Iyz = -0.22 kg.m2 I = [13.31 0.37 0.39 0.37 14.22 -0.22 0.39 -0.22 11.2]; % Momentos principais de inercia [v, e] = eig(I); I1 = e(1,1); I2 = e(2,2); I3 = e(3,3); % inercia da roda Iw = 12/(45000*(pi/180)) ; hs = 12; %saturacao ws = 7500*6 ;%7500rpm*6graus/s a2rpm = 1/6; %--------------------------------------------------------------------% Dinamica linearizada: satelite + rodas %---------------------------------------------------------------------% Matriz A A =[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;... (-4*n^ 2*(I2-I3)-n*H o)/I1 0 0 0 0 (n*(-I2+I3+I1)-H o)/I1;... 0 3*n^ 2*(-I1+I3)/I2 0 0 0 0;... 0 0 (n^ 2*(-I2+I1)-n*H o)/I3 (-n*(I1-I2+I3)+H o)/I3 0 0]; % Matriz B B=[0 0 0 0 0 0 0 0 0 -1/I1 0 0 0 -1/I2 0 0 0 -1/I3]; %Matriz C e D C=eye(6); D=zeros(6,3); 194 Tabela D.2- Programa MATLAB para o projeto LQR do modelo satélite equipado com rodas (cont.). % Calculo da matriz de controlabilidade co = ctrb(A,B); posto co = rank(co); posto A = rank(A); %------------------------------------------------% Lei LQ Regulador %-----------------------------------------------% Matrizes R e Q (ganho) ref. Overby, 2004 rc = 1/((5e-3)^ (2)); % u = 5mNm qc1 = 1/((10*(pi/180))^ (2)); % x = 10 graus qc2 = 1/((1*(pi/180))^ (2)); % xdot = 1 grau/s Rc = rc*diag([1,1,1]); %peso do controle Qc = diag([qc1,qc1,qc1,qc2,qc2,qc2]); %peso do estado % Matriz Kc: ganho do controlador [Kc, S, E] = LQR(A, B, Qc, Rc); %----------------------------------------------% Tracking r(t) %----------------------------------------------% u = Kc*x + v % v = -inv(R)B’s % ds/dt = -[A’-K*B*inv(R)*B’]s + Q*r(t) % r(t) vector reference At = - [A’-S*B*inv(Rc)*B’]; Bt = Qc; Ct = C; Dt = zeros(size(A)); % simulink run(’SimLQEquarsRodasSat’) % end 195 Tabela D.3- Programa MATLAB para o projeto LQR do modelo satélite equipado com bobinas. % Por Gilberto Arantes INPE 2004 % % sintaxe [K] = LQcontrolNcube(B) % Output % K = control matrix gain % Input % B = magnetic field % -----------------------------% LQ Initialization %------------------------------function [K] = LQcontrolNcube(B) % Initial Values Ix = 18.941; Iy = 21.726; Iz = 16.983; n = 1.083e-3; kx = (Iy - Iz)/Ix; ky = (Ix - Iz)/Iy; kz = (Iy - Ix)/Iz; % The Goeomagnetic field Bx o = B(1); By o = B(2); Bz o = B(3); A = [ 0 0 0 1 0 0; 0 0 0 0 1 0; 0 0 0 0 0 1; -4*kx*n^ 2 0 0 0 0 (1-kx)*n; 0 -3*ky*n^ 2 0 0 0 0; 0 0 -kz*n^ 2 -(1-kz)*n^ 2 0 0]; % The input matrix for the linearied system B = [zeros(3,3); 0 Bz o/(2*Ix) -By o/(2*Ix); -Bz o/(2*Iy) 0 Bx o/(2*Iy); By o/(2*Iz) -Bx o/(2*Iz) 0]; % LQ weighting matrices Q = diag([1 1 1 0 0 0])*inv(10*pi/180)^ 2; P = diag([1 1 1])*inv(.1)^ 2; K = -lqr(A,B,Q,P); % Coil Parameters Nx = 100; Ny = 100; Nz = 100; % nunber of coil windings Ax = 0.075; Ay = 0.075; Az = 0.075; % coil area [m^ 2] Rx = 20; Ry = 20; Rz = 20; % Coil resistance [Ohm] M = Nx*Ax; %end 196 Tabela D.4- Programa MATLAB para o projeto LQG do modelo satélite equipado com rodas. % Por Gilberto Arantes @INPE 2003 % % Gera as matrizes para o projeto do sistema de controle LQG % Calculo da equacao matricial de Ricatti e implementacao do LQG % Matrizes A, B, C e D do sistema d/dt(x) = Ax + Bu % y = Cx +du % Inicializacao clc clear all close all deg2rad = pi/180; rad2deg = 180/pi; mi=3.986e14; % moviemnto medio em (m) r=(6378+750)*1000; % raio da orbita circular em (m) n=sqrt(mi/(r^ 3)); % movimento orbital medio em (rad/s) H o = n*0.01; % Momento angular nominal kgm2/s % Momentos e produtos de inercia (14/10/2003 de Sebastiao) % Ixx = 13.31 kg.m2 Ixy = 0.37 kg.m2 % Iyy = 14.22 kg.m2 Ixz = 0.39 kg.m2 % Izz = 11.2 kg.m2 Iyz = -0.22 kg.m2 I = [13.31 0.37 0.39 0.37 14.22 -0.22 0.39 -0.22 11.2]; % Momentos principais de inercia [v,e] = eig(I); I1 = e(1,1); I2 = e(2,2); I3 = e(3,3); % inercia da roda Iw = 12/(45000*(pi/180)) ; hs = 12; %saturacao ws = 7500*6 ;%7500rpm*6graus/s a2rpm = 1/6; %-------------------------------------------------------% Dinamica linearizada: satelite + rodas %-------------------------------------------------------% matriz A A = [0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;... (-4*n^ 2*(I2-I3)-n*H o)/I1 0 0 0 0 (n*(-I2+I3+I1)-H o)/I1;... 0 3*n^ 2*(-I1+I3)/I2 0 0 0 0;... 0 0 (n^ 2*(-I2+I1)-n*H o)/I3 (-n*(I1-I2+I3)+H o)/I3 0 0]; % Matriz B B=[0 0 0 0 0 0 0 0 0 -1/I1 0 0 0 -1/I2 0 0 0 -1/I3]; %Matriz C e D C=eye(6); D=zeros(6,3); 197 Tabela D.5- Programa MATLAB para o projeto LQG do modelo satélite equipado com rodas (cont.). % Calculo da matriz de controlabilidade co = ctrb(A,B); posto co = rank(co); posto A = rank(A); %--------------------------------------------------% Lei LQ Regulador %--------------------------------------------------Matria R e Q (ganho) rc = 1/((5e-3)^ (2)); % u = 5mNm qc1 = 1/((10*(pi/180))^ (2)); % x = 10 graus qc2 = 1/((1*(pi/180))^ (2)); % xdot = 1 grau/s Rc = rc*diag([1,1,1]); %peso do controle Qc = diag([qc1,qc1,qc1,qc2,qc2,qc2]); %peso do estado % Matriz Kc: ganho do controlador [Kc,S,E] = LQR(A,B,Qc,Rc); %--------------------------------------------------% Tracking r(t) %--------------------------------------------------% u = Kc*x + v % v = -inv(R)B’s % ds/dt = -[A’-K*B*inv(R)*B’]s + Q*r(t) % r(t) vector reference At = - [A’-S*B*inv(Rc)*B’]; Bt = Qc; Ct = C; Dt =zeros(size(A)); %-------------------------------------------------%--------------------------------------------------% Implementation Kalman filter % dynamic d/dt(x) = Ax + Bu + Gw % y = Cx + du + v % state estimation % d/dt(xe) = Axe + Bu + Kf(delta - Cxe - Du) delta = yn-ym % d/dt(xe) = Axe + B(-kcxe) + Kf(delta - Cxe - Du) % covariance noise sensor: Evv’ = Rk % covariance noise dynamic Eww’ = Qk G = eye(6,6); sigma angle = (0.5*deg2rad)^ 2; % desvio padrao dos sensor de orientacao sigma rate = (0.05*deg2rad)^ 2; % desvio padrao dos giros Qk = eye(6)*(15*deg2rad)^ 2; % angle +/-(0.5 graus ) Rk = [eye(3)*sigma angle zeros(3); zeros(3) eye(3)*sigma rate]; [Kf,P,E2] = lqe(A,G,C,Qk,Rk); % Structure of compensator % Ac = A - B*Kc - Kf*C; % Bc = Kf; % Simulink run(’SimLQGEquarsRodas’) %end 198 APÊNDICE E TOOLBOX ATITUDE Durante o desenvolvimento do trabalho foi desenvolvido um toolbox para atitude, no ambiente MATLAB e SIMULINK, constituindo uma importante contribuição desse trabalho. Esse apêndice tem por objetivo divulgar o uso do toolbox, apresentando alguns exemplos. A versão final do toolbox será disponibilizado para uso do INPE e instituições interessadas em pesquisas envolvendo dinâmica e controle de atitude. A documetação do toolbox, estará disponı́vel na bibiloteca do INPE. O toolbox é divido nos blocos: • Equações do movimento: Cinemática e dinâmica • Órbita • Ambiente • Controle • Visualização O bloco Equações do movimento contêm as equações da cinemática em diferentes representações: ângulos de Euler e quaternions. O modelo dinâmico da atitude (corpo rı́gido) e o modelo dinâmico do gyrostat (extensão das equações do corpo rı́gido envolvendo movimento interno de rodas). O bloco Órbita contém modelo de órbita Kepleriana, cálculo da latitude e longitude sub-satélite e altura do veı́culo e transformações entre os referenciais envolvidos no estudo de atitude e órbita. O bloco Ambiente contém o modelo do torque de gradiente de gravidade, modelo do campo magnético IGRF, e modelo do torque devido ao momento dipolo residual do satélite. O bloco Controle contém os controladores desenvolvidos nesse trabalho: LQR, LQR tracking, LQG, PID e controladores baseados em energia (Energy based control ). 199 O bloco Visualização contém a visualização/animação da atitude do veı́culo. E.1 E.1.1 Exemplos Equações do Movimento A Figura (E.1) mostra um dos blocos que fazem parte do bloco equações do movimento, o bloco mostrado integra as equações do movimento de um corpo rı́gido contendo rodas (gyrostat), modelo completo, não linear. No ambiente SIMULINK existem vários integradores disponı́veis de passo fixo e/ou variável. A janela mostra os parâmetros de entrada. (a) Modelo do Gyrostat (b) Parâmetros de entrada Figura E.1- Bloco Gyrostat. E.1.2 Ambiente O bloco mostrado na Figura (E.2) calcula o modelo completo, não linear, do torque de gradiente de gravidade. A janela mostra os parâmetros de entrada. 200 (a) Modelo do gradiente de gravidade (b) Parâmetros de entrada Figura E.2- Bloco gradiente de gravidade. 201 Livros Grátis ( http://www.livrosgratis.com.br ) Milhares de Livros para Download: Baixar livros de Administração Baixar livros de Agronomia Baixar livros de Arquitetura Baixar livros de Artes Baixar livros de Astronomia Baixar livros de Biologia Geral Baixar livros de Ciência da Computação Baixar livros de Ciência da Informação Baixar livros de Ciência Política Baixar livros de Ciências da Saúde Baixar livros de Comunicação Baixar livros do Conselho Nacional de Educação - CNE Baixar livros de Defesa civil Baixar livros de Direito Baixar livros de Direitos humanos Baixar livros de Economia Baixar livros de Economia Doméstica Baixar livros de Educação Baixar livros de Educação - Trânsito Baixar livros de Educação Física Baixar livros de Engenharia Aeroespacial Baixar livros de Farmácia Baixar livros de Filosofia Baixar livros de Física Baixar livros de Geociências Baixar livros de Geografia Baixar livros de História Baixar livros de Línguas Baixar livros de Literatura Baixar livros de Literatura de Cordel Baixar livros de Literatura Infantil Baixar livros de Matemática Baixar livros de Medicina Baixar livros de Medicina Veterinária Baixar livros de Meio Ambiente Baixar livros de Meteorologia Baixar Monografias e TCC Baixar livros Multidisciplinar Baixar livros de Música Baixar livros de Psicologia Baixar livros de Química Baixar livros de Saúde Coletiva Baixar livros de Serviço Social Baixar livros de Sociologia Baixar livros de Teologia Baixar livros de Trabalho Baixar livros de Turismo