UNIVERSIDADE FEDERAL DO PARÁ
INSTITUTO DE TECNOLOGIA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA ELÉTRICA
ANTONIO DA SILVA SILVEIRA
TÉCNICAS DE CONTROLE LQG/LTR E FUZZY APLICADAS AO VEÍCULO AÉREO
NÃO TRIPULADO: AEROSONDE
DM: 02 / 2008
UFPA / ITEC / PPGEE
Campus Universitário do Guamá
Belém-Pará-Brasil
2008
UNIVERSIDADE FEDERAL DO PARÁ
INSTITUTO DE TECNOLOGIA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA ELÉTRICA
ANTONIO DA SILVA SILVEIRA
TÉCNICAS DE CONTROLE LQG/LTR E FUZZY APLICADAS AO VEÍCULO AÉREO
NÃO TRIPULADO: AEROSONDE
Dissertação submetida à
Banca
Examinadora
do
Programa de Pós-Graduação
em Engenharia Elétrica da
UFPA para a obtenção do
Grau
de
Mestre
em
Engenharia Elétrica
UFPA / ITEC / PPGEE
Campus Universitário do Guamá
Belém-Pará-Brasil
2008
________________________________________________________________________
S587t
Silveira, Antonio da Silva
Técnicas de controle LQG/LTR e fuzzy aplicadas ao veículo aéreo não
tripulado : Aerosonde / Antonio da Silva Silveira; orientador, José Augusto
Lima Barreiros.- 2008.
Dissertação (Mestrado) – Universidade Federal do Pará, Instituto de
Tecnologia, Programa de Pós-Graduação em Engenharia Elétrica, Belém,
2008.
1. Controle automático. 2. Aviões - sistemas de controle 3. Sistemas difusos.
I. Título.
CDD – 22. ed. 629.89
UNIVERSIDADE FEDERAL DO PARÁ
INSTITUTO DE TECNOLOGIA
PROGRAMA DE PÓS-GRADUAÇÃO EM ENGENHARIA ELÉTRICA
TÉCNICAS DE CONTROLE LQG/LTR E FUZZY APLICADAS AO VEÍCULO AÉREO NÃO
TRIPULADO: AEROSONDE
AUTOR: ANTONIO DA SILVA SILVEIRA
DISSERTAÇÃO DE MESTRADO SUBMETIDA À AVALIAÇÃO DA BANCA EXAMINADORA
APROVADA PELO COLEGIADO DO PROGRAMA DE PÓS-GRADUAÇÃO EM
ENGENHARIA ELÉTRICA DA UNIVERSIDADE FEDERAL DO PARÁ E JULGADA
ADEQUADA PARA OBTENÇÃO DO GRAU DE MESTRE EM ENGENHARIA ELÉTRICA NA
ÁREA DE SISTEMAS DE ENERGIA
APROVADA EM: 25 / 02 / 2008
BANCA EXAMINADORA:
_________________________________________________________________________________
Prof. Dr. JOSÉ AUGUSTO LIMA BARREIROS
(ORIENTADOR – UFPA)
_________________________________________________________________________________
Prof. Dr. CARLOS TAVARES DA COSTA JÚNIOR
(CO-ORIENTADOR – UFPA)
_________________________________________________________________________________
Prof. Dr. ANDRÉ MAURÍCIO DAMASCENO FERREIRA
(MEMBRO – CEFET / PA)
_________________________________________________________________________________
Profa. Dra. BRIGIDA RAMATI PEREIRA DA ROCHA
(MEMBRO – UFPA)
_________________________________________________________________________________
Prof. Dr. JORGE ROBERTO BRITO DE SOUZA
(MEMBRO – UFPA)
_________________________________________________________________________________
Prof. Dr. ORLANDO FONSECA SILVA
(MEMBRO – UFPA)
_________________________________________________________________________________
Prof. Dr. WALTER BARRA JÚNIOR
(MEMBRO – UFPA)
VISTO:
_________________________________________________________________________________
Prof. Dr. EVALDO GONÇALVES PELAES
(COORDENADOR DO PPGEE / ITEC / UFPA)
UFPA / ITEC / PPGEE
Dedicatória
Dedico este trabalho ao meu avô, Euzébio Orlando da Mota Silveira, que
infelizmente se foi quando eu ainda estava na graduação. Tenho certeza que, onde
quer que ele esteja, sem dúvida, deve estar muito contente pela conclusão deste
trabalho... Tanto quanto eu fico quando me lembro dele.
Agradecimentos
Agradeço a dedicação e apoio dado por minha esposa, Aline Silveira, que
contribuiu com este trabalho, desde sua concepção, quando iniciamos o mesmo
durante as disciplinas do curso de mestrado. Trabalhamos juntos enquanto a idéia
ainda estava se consolidando, e chegamos a elaborar juntos, a primeira versão, na
forma de um artigo, que posteriormente se tornou esta dissertação.
Agradeço, também, aos meus pais, Antonio e Leda, e a minha irmã, Luciana,
pela motivação, incentivo, patrocínio e muitas vezes, até auxílio técnico e teórico,
quando eram obrigados a ler e ouvir durante horas, minhas histórias de dissertação
e aviação, e me ajudar a conduzir alguma tarefa relacionada ao meu trabalho.
Todos os meus amigos, de alguma forma, me apoiaram e incentivaram,
mesmo quando nem se quer tinham ciência disso, mas em especial, gostaria de
agradecer aos amigos Cleison Daniel Silva, Rafael Bayma e ao professor e amigo,
Max Rothe-Neves, que muito me ensinaram e participaram na minha formação.
Agradeço também, aos amigos aeromodelistas do Aeroclube do Pará, com quem
aprendi bastante durante os divertidos finais de semana em que brincamos de piloto,
mesmo que remoto. Além deles, agradecimento especial ao professor e orientador
Barreiros, meu co-orientador Tavares, Prof. Orlando (Nick), Prof. Jorge Brito e ao
Prof. Walter Barra – eles são os principais responsáveis pela minha formação na
área de sistemas e controle.
Finalizando, gostaria de agradecer o patrocínio da CAPES, que financiou
meus estudos durante o curso do mestrado, e juntamente, agradecer a equipe do
Programa de Pós-Graduação da Universidade Federal do Pará.
Sumário
1 Introdução...............................................................................................................1
2 Física do processo, geodésia e navegação .........................................................5
2.1 Cinemática e dinâmica de corpos rígidos ..........................................................5
2.1.1 Cinemática vetorial......................................................................................5
2.1.2 Dinâmica de corpos rígidos.......................................................................12
2.2 Geodésia .........................................................................................................16
2.3 Navegação terrestre ........................................................................................19
2.3.1 Frames de referência e sistemas de coordenadas....................................19
2.3.2 Sistema sexagesimal e posição na Terra..................................................20
2.3.3 Método de navegação adotado.................................................................21
2.4 Aerodinâmica básica........................................................................................22
2.4.1 Aerodinâmica da seção de um aerofólio ...................................................22
3 Aerosim 1.2, introdução aos UAVs e ao Aerosonde .........................................26
3.1 Aerosim 1.2......................................................................................................26
3.1.1 Interface ....................................................................................................27
3.2 Introdução aos UAVs .......................................................................................29
3.3 Aerosonde UAV ...............................................................................................31
3.3.1 Estudo da dinâmica do Aerosonde ...........................................................32
3.3.1.1 Dinâmica longitudinal..........................................................................38
3.3.1.2 Dinâmica lateral / direcional................................................................43
4 Visão geral dos sistemas de controle de aeronaves ........................................48
4.1 Especificações de pólos/modos.......................................................................49
4.2 Sistemas de aumento de estabilidade (SAS)...................................................50
4.2.1 Pitch SAS ..................................................................................................51
4.2.2 SAS Lateral-Direcional ou Yaw Damper ...................................................52
4.3 Sistemas de aumento de controle (CAS).........................................................53
4.3.1 Pitch-Rate CAS .........................................................................................54
4.4 Pilotos Automáticos .........................................................................................55
4.4.1 Pitch-Attitude Hold.....................................................................................56
4.4.2 Altitude Hold ..............................................................................................56
4.4.3 Speed Hold ...............................................................................................57
4.4.4 Roll Angle Hold..........................................................................................58
4.5 Técnicas de controle robusto multivariável ......................................................58
4.5.1 Robustez e estabilidade ............................................................................59
4.5.2 Especificações de desempenho no domínio da freqüência ......................60
4.5.3 LQG...........................................................................................................61
4.5.3.1 LQR ....................................................................................................63
4.5.3.2 Filtro de Kalman..................................................................................64
4.5.4 LQG/LTR...................................................................................................67
4.6 Sistema de Controle Lateral do Aerosonde UAV.............................................69
4.6.1 LQG/LTR Lateral do Aerosonde................................................................73
4.7 Sistema de Controle Longitudinal do Aerosonde UAV ....................................80
4.7.1 LQG/LTR Longitudinal do Aerosonde .......................................................82
5 Controladores fuzzy de Altitude, de Direção e o Sistema de Navegação .......89
5.1 Controlador fuzzy de Altitude...........................................................................89
5.2 Controlador fuzzy de direção/heading .............................................................92
5.3 Sistema de navegação ....................................................................................94
5.4 Simulações com o piloto automático completo ................................................96
5.4.1 Simulação 1: na região de estabilidade e sem perturbações ....................97
5.4.1.1 Planejamento da missão ..................................................................101
5.4.1.2 Teste dos controladores com o modelo não-linear ...........................105
5.4.1.3 Simulação da rota: Sobrevôo em Belém ..........................................109
5.4.2 Simulação 2: com variação de parâmetros, perturbações e fora da região
de estabilidade programada.............................................................................112
6 Conclusão Geral e Sugestões para novos trabalhos......................................116
Referências Bibliográficas ...................................................................................119
APÊNDICE A - Instalação do Aerosim 1.2 para MATLAB 6 ou superior...........121
APÊNDICE B – Controlador LQG/LTR Lateral do Aerosonde ...........................123
APÊNDICE C – Controlador LQG/LTR Longitudinal do Aerosonde..................127
APÊNDICE D – Controlador fuzzy de Altitude ....................................................132
APÊNDICE E – Controlador fuzzy de Direção.....................................................134
APÊNDICE F – Gerador de waypoints .................................................................136
APÊNDICE G – Sistema de Navegação ...............................................................137
ANEXO A – Software de Configuração do Aerosonde UAV ..............................139
ANEXO B – Software de Trim e Linearização do Aerosonde ............................145
Lista de Figuras
2. 1: Derivada de um vetor em um frame em rotação. ................................................7
2. 2: Velocidade e aceleração com frames em movimento. ........................................8
2. 3: Geoid e definições de altura..............................................................................18
2. 4: Modelo esferoidal oblato da Terra.....................................................................18
2. 5: Longitude e latitude em relação ao Equador e o meridiano de Greenwich. ......21
2. 6: Seção de um aerofólio e principais parâmetros. ...............................................23
2. 7: Definição das superfícies de controle, de ângulos e eixos aerodinâmicos........24
3. 1: Interface da biblioteca de blocos do Aerosim 1.2 no Simulink...........................27
3. 2: Fotografia e modelo gráfico tridimensional do Aerosonde UAV. .......................31
3. 3: Diagrama de simulação no Simulink do Aerosonde tipo caixa preta.................32
3. 4: Resposta temporal do ângulo de bank. .............................................................33
3. 5: Resposta temporal de heading..........................................................................33
3. 6: Resposta temporal de φ para entrada nos ailerons de -0,0073 rad..................34
3. 7: Resposta temporal de φ do modelo não-linear com os valores de entrada
trimados. ............................................................................................................36
3. 8: Resposta temporal da velocidade relativa (airspeed) para o ponto de operação
em 23 m/s. .........................................................................................................37
3. 9: Resposta temporal do ângulo θ para o ponto de operação de trimagem. ........37
3. 10: Resposta temporal de α e θ ao doublet no elevator e modo natural shortperiod estimulado. ..............................................................................................39
3. 11: Resposta temporal de h e Va ao doublet no elevator......................................40
3. 12: Resposta temporal de α e θ ao doublet no throttle e oscilação em θ devido
ao phugoid. ........................................................................................................41
3. 13: Resposta temporal de h e Va ao doublet no throttle. Phugoid estimulado. .....42
3. 14: Resposta de β e φ ao doublet aplicado no rudder. .......................................44
3. 15: Resposta de β e ψ ao doublet aplicado no rudder........................................45
3. 16: Resposta de φ e ψ ao doublet aplicado no rudder. .......................................45
4. 1: Estrutura de um Pitch SAS com realimentação do ângulo de ataque e pitch rate.
...........................................................................................................................51
4. 2: Estrutura de um Yaw Damper. ..........................................................................52
4. 3: Perda de sustentação durante o roll influenciando na dinâmica longitudinal e
lateral. ................................................................................................................54
4. 4: Diagrama de um possível pitch-rate CAS..........................................................55
4. 5: Diagrama de blocos de um piloto automático tipo Pitch-Attitude Hold. .............56
4. 6: Diagrama de blocos de um piloto automático de altitude com sub-controlador de
ângulo de pitch...................................................................................................57
4. 7: Diagrama de blocos de um Speed Hold e Altitude Hold na mesma malha de
controle. .............................................................................................................57
4. 8: Estrutura do controlador LQG. ..........................................................................62
4. 9: Estrutura da planta com observador de estados. ..............................................66
4. 10: Diagrama de blocos do controlador LQG/LTR (Zarei et al.; 2006). .................68
4. 11: Ganhos principais do sistema original da dinâmica lateral do Aerosonde.......72
4. 12: Ganhos principais balanceados em DC e com inclinação de -20 dB/década..72
4. 13: Magnitude de
1
m (ω )
, para limitar σ nas altas freqüências................................73
4. 14: Ganhos principais de
CΦ ( s) L
com Filtro de Kalman e rf = 10 2 . .........................74
4. 15: Resposta ao degrau aplicado em rφ do sistema em M.F. com Filtro de Kalman.
...........................................................................................................................75
4. 16: Erro de estimação dos estados tendendo a zero, exceto por ψ na linha preta
tracejada. ...........................................................................................................75
4. 17: Diagrama de simulação do Simulink e análise do erro de estimação de
estados...............................................................................................................76
4. 18: Ganhos principais do sistema com LQG/LTR. ................................................78
4. 19: Resposta ao degrau do sistema com LQG/LTR e sinal de controle................78
4. 20: Ganhos principais inferior e superior do sistema longitudinal original. ............81
4. 21: Ganhos principais do sistema aumentado com pré-compensador. .................81
4. 22: Ganhos principais do sistema aumentado em M.A. com o Filtro de Kalman...83
4. 23: Resposta temporal ao degrau aplicado em rθ do sistema em M.F. com o Filtro
de Kalman. .........................................................................................................83
4. 24: Resposta temporal ao degrau aplicado em rV a do sistema em M.F. com o Filtro
de Kalman. .........................................................................................................84
4. 25: Erros de estimação dos estados do sistema com o Filtro de Kalman. ............84
4. 26: Ganhos principais do sistema com LQG/LTR. ................................................86
4. 27: Resposta ao degrau aplicado em rθ do sistema em M.F. com LQG/LTR e o
sinal de controle. ................................................................................................87
4. 28: Resposta ao degrau aplicado em rVa do sistema em M.F. com LQG/LTR e o
sinal de controle. ................................................................................................87
5. 1: Diagrama de blocos do controlador fuzzy de Mamdani.....................................90
5. 2: Funções de pertinência de entrada do controlador. ..........................................91
5. 3: Funções de pertinência de saída do controlador...............................................91
5. 4: Funções de pertinência de entrada do controlador. ..........................................93
5. 5: Funções de pertinência da saída do controlador...............................................93
5. 6: Representação gráfica do método de triangulação para obtenção da trajetória
ψ ref .....................................................................................................................95
5. 7: Diagrama completo de simulação no Simulink..................................................99
5. 8: Diagrama interno do bloco 1. ..........................................................................100
5. 9: Diagrama interno do bloco 2. ..........................................................................100
5. 10: Diagrama interno do bloco 3. ........................................................................100
5. 11: Diagrama interno do bloco 4. ........................................................................100
5. 12: Diagrama interno do bloco 5. ........................................................................101
5. 13: Diagrama interno do bloco 7. ........................................................................101
5. 14: Rota que deverá ser percorrida passando por cinco waypoints. ...................102
5. 15: Rota gerada pelo waypoint_gen.m. ...............................................................104
5. 16: Resposta temporal do controle de θ e φ com zoom nos primeiros 4 s de
simulação. ........................................................................................................105
5. 17: Resposta temporal de Va para o teste dos controladores com asas niveladas.
.........................................................................................................................105
5. 18: Resposta temporal da simulação com doublet aplicado em φref . ..................106
5. 19: Resposta temporal de h e intervalo
marcando o momento em que θ reduz
e Va aumenta. .................................................................................................108
5. 20: Resposta temporal de φ, θ , β , e α no controle de altitude e direção. ..........108
5. 21: Resposta temporal de Va e ψ no controle de altitude e direção. .................109
5. 22: Gráfico 3-D da rota percorrida pelo Aerosonde UAV no Sobrevôo em Belém.
.........................................................................................................................110
5. 23: Gráfico 2-D da rota percorrida pelo Aerosonde UAV no Sobrevôo em Belém.
.........................................................................................................................110
5. 24: Variação da altitude ao longo da simulação. .................................................111
5. 25: Variação de ψ ao longo da simulação. .........................................................111
5. 26: Variação de Va ao longo da simulação.........................................................111
5. 27: Velocidades dos ventos sul e oeste e da turbulência usada na simulação. ..112
5. 28: Rota 3-D percorrida pelo Aerosonde submetido a perturbações atmosféricas e
variação de parâmetros. A altitude dos waypoints entre 250 e 300 m . ...............113
5. 29: Gráfico 2-D da Simulação 1 e 2 juntas para comparação entre simulação com
e sem perturbações e variações de parâmetros. .............................................114
5. 30: Direção e Altitude foram afetadas pelos ventos e turbulências. ....................114
5. 31: Va e Vamédia ao longo da simulação com perturbações e variação de parâmetros.
.........................................................................................................................115
Lista de Tabelas
2. 1: Convenções de sinais para as superfícies de controle aerodinâmicas. ............25
3. 1: Pólos, fatores de amortecimento, freqüências e constantes de tempo do modelo
longitudinal. ........................................................................................................42
3. 2: Pólos, fatores de amortecimento, freqüências e constantes de tempo da
dinâmica lateral. .................................................................................................47
4. 1: Tipos comuns de sistemas de aumento de estabilidade, de controle e pilotos
automáticos........................................................................................................49
4. 2: Pólos de malha fechada do sistema com o Filtro de Kalman ............................76
4. 3: Em cinza, os pólos M.F. do Aerosonde aumentados pelo pré-compensador e
modificados pelo controlador. Em branco, os pólos adicionados pelo LQG/LTR.
...........................................................................................................................79
4. 4: Pólos de M.F. do sistema aumentado com o Filtro de Kalman. ........................85
4. 5: Em cinza, os pólos M.F. do Aerosonde aumentados pelo pré-compensador e
modificados pelo controlador. Em branco, os pólos adicionados pelo LQG/LTR.
...........................................................................................................................88
Resumo
Este trabalho é um estudo sobre sistemas de controle de vôo de aviões, que aborda
desde um resumo teórico sobre a física do processo, aerodinâmica, navegação
terrestre e geodésia, até o projeto de controladores LQG/LTR e fuzzy, para compor
um piloto automático testado em simulações com um modelo não-linear da planta, o
veículo aéreo não tripulado, Aerosonde UAV. Além disto, um sistema de navegação
é projetado para simular o planejamento e seguimento de uma rota de vôo préestabelecida em duas simulações, uma com perturbações meteorológicas e variação
de parâmetros do modelo, e outra em condições ideais. O que garante o bom
desempenho deste piloto automático, mesmo quando testado em condições
adversas, tal como em pontos de operação diferentes do qual foi projetado, é a
utilização do método LQG/LTR nos mais baixos níveis da malha de controle,
compondo os controladores longitudinal e lateral do Aerosonde, que atuam
diretamente no comando das principais superfícies de controle aerodinâmicas e na
potência do motor da aeronave. O método LQG/LTR usufrui das propriedades do
Regulador Linear Ótimo Quadrático (LQR) com realimentação total dos estados, que
por sua vez, é garantidamente robusto e estável. Além disto, o método LQG/LTR
também utiliza o Filtro de Kalman como estimador de estados e para filtrar dinâmicas
indesejáveis de alta freqüência, sendo o projeto baseado em análises no domínio da
freqüência dos ganhos principais do sistema. Já na camada intermediária da malha
de controle, são utilizados controladores fuzzy do tipo Mamdani para controlar a
altitude e direção da aeronave através de um mapeamento estático das variáveis
erro de altitude e direção. As saídas destes controladores fuzzy servem de
referência para os controladores LQG/LTR longitudinal e lateral. Na camada mais
alta da malha de controle, está presente o sistema de navegação e coordenação de
vôo, que monitora e determina a trajetória a ser percorrida pela aeronave, passando
estas informações como referência aos controladores fuzzy. Desta forma, o sistema
completo usufrui da fusão de técnicas modernas e inteligentes de controle
procurando atender uma vasta região de operação sem a necessidade de utilizar
ganhos programáveis.
Palavras-chave: controle de vôo, aviões, método LQG/LTR, controle fuzzy, piloto
automático, Aerosonde, UAV, sistema de navegação.
Abstract
This work is a study about aircraft flight control systems that covers from a theoretical
resume about the process physics to the design of LQG/LTR and fuzzy controllers, to
build an autopilot tested in simulations with a non-linear model of the plant, the
unmanned aerial vehicle, Aerosonde UAV. Furthermore, a navigation system is
designed for the planning and tracking of a pre-established flight route in two
simulations, one with meteorological disturbances and parameter variations of the
model, and another in ideal conditions. What guarantees a good performance of this
autopilot, even when tested in adverse conditions, as in different operation points
from where it was designed for, is the use of the LQG/LTR method inside of low
levels of the control loop, composing the longitudinal and lateral Aerosonde
controllers, that acts directly over the principal aerodynamic command surfaces and
the engine power of the aircraft. The LQG/LTR method usufructs of properties from
the Linear Quadratic Regulator (LQR) with full state feedback, with it´s guaranteed
robustness and stability. Along with that, the LQG/LTR method also uses the Kalman
Filter as a state estimator and to filter high frequency undesirable dynamics, being
the design process based on frequency domain analisys of the principal gains of the
system. And inside of an intermediary level of the control loop, fuzzy controllers of
Mamdani kind are used to control altitude and direction of the aircraft, by means of a
stactic mapping of altitude and heading error variables. The outputs of these fuzzy
controllers are used as reference inputs for the LQG/LTR longitudinal and lateral
controllers. In the higher level of the control loop, is where the navigation and flight
coordination system are present, monitoring and determining the right trajactory to be
tracked by the aircraft, passing that information as reference inputs for the fuzzy
controllers. In this way, the complete system usufructs of the fusion between modern
and intelligent control techniques, looking forward to attend a wide operational range
without the need to use scheduled gains.
Keywords: flight control, aircraft, LQG/LTR method, fuzzy control, autopilot,
Aerosonde, UAV, navigation system.
Capítulo 1
Introdução
Desde o início da história da aviação, deu-se ênfase em construir aviões
controláveis pela pilotagem mais do que pela estabilidade inerente dos mesmos. As
dificuldades em controlar os primeiros aviões assim como o crescimento na duração
e distâncias dos vôos, rapidamente, motivaram o desenvolvimento de sistemas de
controle automático adequados ao ramo da aviação. Atualmente, tanto o tempo, a
distância e as fases de vôo são as mais diversas possíveis, exigindo constante
adaptação dos sistemas de controle da aeronave.
As técnicas e aplicações da engenharia de controle e sistemas na solução de
problemas no ramo da aviação são bem diversas, podendo-se encontrar na literatura
especializada métodos baseados em teoria de controle clássico, moderno,
inteligente e pela fusão dessas teorias.
Técnicas de controle clássico, sozinhas, apresentam certas limitações no
controle de aeronaves. Isto se deve ao grande número de variáveis de entrada e
saída (muitas vezes fortemente acopladas), à necessidade de boa experiência e
intuição para definir a estrutura dos controladores, ao fechamento ordenado e
sucessivo das malhas de realimentação e à sintonia de boa parte dos controladores
por tentativa e erro.
A teoria de controle moderno tem tido impacto significativo na indústria
aeronáutica nos últimos anos, utilizando técnicas como: alocação de pólos, modelo
de referência, inversão dinâmica, LQG/LTR (Linear Quadratic Gaussian/LoopTransfer-Recovery) e LQ com realimentação da saída (Stevens e Lewis; 2003). Uma
das principais vantagens em utilizar tais técnicas é que o projeto é baseado
diretamente no modelo por variáveis de estado, que normalmente contém mais
informação sobre o sistema do que a descrição por função de transferência baseada
apenas na relação entrada-saída (modelo tipo caixa preta). Outra vantagem é a
descrição do sistema por equações matriciais, que ao contrário das técnicas
clássicas, podem computar todos os ganhos de controle simultaneamente e todas as
malhas de realimentação podem ser fechadas ao mesmo tempo, permitindo projetos
mais rápidos e diretos.
1
Técnicas de controle inteligente aplicando, por exemplo, lógica fuzzy, também
têm sido usadas em trabalhos recentes, tais como: Sugeno et al. (1993), no controle
de um helicóptero não tripulado; Cavalcante (1994), na elaboração de um sistema de
navegação para helicóptero não tripulado; Doitsidis et al. (2004), Mileva at al. (2005)
e Silveira et al. (2007), com simulações de pequenos aviões tripulados e não
tripulados utilizando controladores fuzzy do tipo Mamdani (Passino e Yurkovich;
1998). Nesses últimos, principalmente, apesar de funcionarem e fornecerem boas
respostas temporais, os controladores são sintonizados por tentativa e erro, ou
baseando-se no funcionamento de controladores clássicos para sintonizar
controladores fuzzy do tipo Mamdani, o que remete aos problemas de definição das
estruturas dos controladores e fechamentos sucessivos de malhas de realimentação,
sem provas matemáticas das margens de estabilidade e de robustez do sistema
completo.
As técnicas mais recentes procuram fundir as teorias clássicas, modernas e
inteligentes na tentativa de usufruir de suas vantagens e suprir as necessidades com
o que a outra possa oferecer. Neste trabalho, o estudo é voltado à fusão das
técnicas inteligentes e modernas utilizando sistemas fuzzy do tipo Mamdani que
analisam as saídas da planta e comandam compensadores LQG/LTR que garantem
a robustez e estabilidade do sistema levando em consideração margens de
incertezas do modelo e ruídos de medição. Outra característica interessante do
LQG/LTR é que este dispensa grande intuição sobre a estrutura do compensador,
fornecendo-a de forma direta.
Na indústria aeronáutica costuma-se utilizar gain-scheduling (ou ganhos
programáveis) devido a variações nas fases de vôo e alterações nos mais diversos
parâmetros da aeronave. Gain-scheduling se baseia em determinar ganhos (ou
parâmetros) do controlador para vários pontos de operação da aeronave e usando
um sistema supervisório, selecionar estes ganhos quando a aeronave estiver no
respectivo ponto de operação. No entanto, dependendo da quantidade significativa
destas variações, o trabalho de obtenção dos modelos lineares locais para o projeto
dos controladores e mais o projeto em si, pode ser exaustivo, necessitando-se em
alguns casos, de centenas de projetos. Com o uso do LQG/LTR pretende-se utilizar
um único modelo linear para os projetos levando em consideração as margens de
incertezas do modelo e assim garantir sua utilização em uma grande faixa de
operação.
2
Para este trabalho, a aeronave utilizada nas simulações não possui um
envelope (ou faixa de operação) velocidade-altitude muito grande. Ela é um micro
veículo aéreo não tripulado (micro UAV, do inglês, Unmanned Aerial Vehicle) que
opera apenas em velocidades subsônicas. Mesmo assim, executa boa parte das
fases de vôo como qualquer outra aeronave tripulada, como: decolagem, subida até
altitude de cruzeiro, cruzeiro, descida, aproximação e pouso. Logo, por gainscheduling simplesmente, a varredura de um grande número de pontos de operação
e uma grande tabela de ganhos teria de ser obtida, diferentemente da solução por
LQG/LTR junto com controle fuzzy que será abordada.
Além da preocupação em cobrir uma vasta região de operação da aeronave,
aviões são sistemas multivariáveis bastante complexos, com um grande número de
variáveis que não são medidas pelos sensores e com forte acoplamento dos canais
de entrada e saída, sendo então necessário estimar algumas variáveis e garantir o
desacoplamento dos canais citados sob condições de restrições que otimizem o
esforço de controle e minimizem os efeitos das dinâmicas indesejadas. O LQG/LTR
utiliza então o Filtro de Kalman para estimar as variáveis e permitir a realimentação
total de estados para o LQR (Regulador Linear Ótimo Quadrático, do inglês, Linear
Quadratic Regulator).
A parte final do trabalho define o piloto automático, interliga os sistemas de
controle e estabilização da aeronave (sistemas de baixo nível com LQG/LTR) aos
controladores fuzzy de Mamdani e aos sistemas de navegação e coordenação de
vôo. O sistema de navegação é responsável por informar dados de posição do
veículo no globo terrestre e o de coordenação, por permitir o planejamento dos vôos
e determinar e passar comandos aos sistemas de baixo nível.
Para simular tais sistemas, como: a atmosfera, campo magnético, forma do
planeta Terra, gravidade, Sistema de Posicionamento Global (GPS, do inglês, Global
Positioning System) e o modelo não-linear do avião com seis graus de liberdade (6DOF, do inglês, Six Degrees Of Freedom), será utilizada uma ferramenta para o
software MATLAB (http://www.mathworks.com) chamada Aerosim Blockset, da
empresa U-Dynamics (http://www.u-dynamics.com).
3
O texto deste trabalho, além deste capítulo introdutório, foi estruturado da
seguinte forma:
Capítulo 2: revisão de conceitos sobre cinemática e dinâmica de aeronaves,
aerodinâmica básica, forças e momentos que agem sobre os aviões, geodésia e
navegação terrestre.
Capítulo 3: apresentação do ambiente de simulação e da ferramenta Aerosim
1.2; introdução aos UAVs e ao Aerosonde UAV que será utilizado nas simulações;
apresentação do modelo matemático linear do Aerosonde UAV, estudo de sua
dinâmica e análise dos modos naturais.
Capítulo 4: visão geral dos sistemas de controle, de estabilização e pilotos
automáticos; revisão teórica de técnicas de controle robusto multivariável e análise
de sistemas multivariáveis no domínio da freqüência; projeto dos controladores
LQG/LTR de controle e estabilização lateral e longitudinal.
Capítulo 5: projeto dos controladores fuzzy de Mamdani de altitude e direção
que comandam os sub-controladores LQG/LTR; desenvolvimento do sistema de
navegação; resultados das simulações com o piloto automático completo;
Capítulo 6: conclusões gerais e sugestões para continuações do trabalho.
4
Capítulo 2
Física do processo, geodésia e navegação
O objetivo deste capítulo é familiarizar o leitor aos tópicos considerados
essenciais ao entendimento da engenharia aeroespacial, dando ênfase a partes
relevantes da mecânica clássica voltada ao tratamento de veículos com seis graus
de liberdade (6-DOF) através de uma representação por ângulos de Euler e
quatérnions. Além desses tópicos, simulações de veículos em movimento sobre e ao
redor da Terra requerem um resumo sobre geodésia, gravitação e navegação.
Finalizando, um estudo sobre aerodinâmica básica fornecerá a base necessária para
compreensão dos capítulos seguintes, sendo considerada indispensável ao
entendimento dos mesmos.
2.1 Cinemática e dinâmica de corpos rígidos
A cinemática e dinâmica de corpos rígidos não considera estruturas flexíveis,
logo, as aeronaves neste trabalho sempre mantém todas as suas partes (não
móveis) na mesma posição relativa. O que, de acordo com Stevens e Lewis (2003),
na maioria dos casos é suficiente para simulações e projetos de sistemas de
controle de vôo quando não se está preocupado em aliviar cargas aerodinâmicas e
controlar modos estruturais.
2.1.1 Cinemática vetorial
A cinemática estuda o movimento de corpos sem levar em conta os
mecanismos causadores do movimento. O movimento de objetos físicos pode ser
descrito por vetores em três dimensões utilizando as seguintes definições:
Frame de Referência: corpo rígido ou conjunto de pontos rigidamente
acoplados usados para estabelecer distâncias e direções, podendo estar em
movimento ou ser um referencial inercial;
Vetor: objeto geométrico abstrato que possui magnitude e direção;
5
Sistema de coordenadas: sistema de medidas para localização de pontos
num espaço definido em um frame de referência.
Notações
Algumas notações importantes também devem ser destacadas, tal como em
Stevens e Lewis (2003):
p A / B ≡ vetor posição do ponto A com respeito ao ponto B
v A / i ≡ velocidade do ponto A em relação ao frame i (Fi )
b
v A / i ≡ derivada de v A / i calculada em Fb
v cA / i ≡ componentes de v A / i no sistema de coordenadas c
b
v cA / i ≡ componentes no sistema de coordenadas c da derivada em Fb
Exemplo dos componentes de um vetor:
⎡ xb ⎤
⎡ vx ⎤
b
⎢
⎥
⎢ ⎥
p A / B = ⎢ yb ⎥ e v = ⎢v y ⎥
⎢⎣ zb ⎥⎦
⎢⎣ vz ⎥⎦
b
são componentes de p e v em um sistema de coordenadas b.
Vetor de velocidade angular
A derivada de um vetor pode ser definida de maneira análoga à definição de
derivada de uma grandeza escalar:
⎡ p (t + δ t ) − p A/ B (t ) ⎤
dp A / B
= lim ⎢ A / B
⎥
δ t →0
dt
δt
⎣
⎦
Se p A / B é um vetor posição, sua derivada é um vetor velocidade somente se a
derivada for calculada no frame no qual B é um ponto fixo.
A Figura 2.1 exemplifica a derivada de um vetor em um frame rotativo com
respeito a um vetor de referência. A derivada de
w tomada no
frame
Fr é não nula se w estiver mudando de direção ou magnitude quando observado de
Fr , sendo independente de qualquer translação entre os dois frames. A mudança de
direção pode ser obtida utilizando o teorema da rotação (Stevens e Lewis; 2003). Na
Figura 2.1, seja ŝ um vetor unitário paralelo ao eixo de rotação no tempo t , para um
6
observador em Fr , w se torna um novo vetor w + δ w no tempo t + δ t devido a
rotação δφ . Calculando-se a rotação com δφ positivo em torno de ŝ , fornece:
δ w ⎛ δφ ⎞
≈ sˆ
×w
δ t ⎜⎝ δ t ⎟⎠
quando δ t → 0 ,
r
( )
= sˆφ × w
w
( )
A quantidade sˆφ compõe o vetor velocidade angular instantânea, ωb / r de Fb , com
respeito ao frame Fr . Se w também muda em magnitude em Fb , este efeito é
considerado na equação:
r
= bw
+ ωb / r × w
w
(2.1)
A equação (2.1) é chamada de equação de Coriolis (Stevens e Lewis; 2003) e é
considerada essencial no desenvolvimento de equações de movimento baseadas
nas leis de Newton, aplicando-se a qualquer quantidade física representada por
vetores sem a necessidade de avaliar as derivadas em função do tempo.
É importante destacar algumas propriedades do vetor velocidade angular: 1) é
um vetor único que relaciona as derivadas de um vetor tomadas em dois frames
diferentes; 2) satisfaz a condição de movimento relativo ωb / a = −ωa / b ; 3) é aditivo
sobre múltiplos frames, ωc / a = ωc / b + ωb / a ; 4) sua derivada é a mesma em ambos os
frames, aω b / a = bω b / a .
Figura 2. 1: Derivada de um vetor em um frame em rotação.
7
Velocidade e Aceleração em frames em movimento
A Figura 2.2 apresenta um ponto P movendo-se com respeito aos frames
Fa e Fb em relação aos respectivos pontos fixos O e Q . Para relacionar as
velocidades e acelerações nos dois frames, primeiro obtêm-se os vetores de posição
e as derivadas em relação ao frame Fa .
rP / O = rQ / O + rP / Q
a
rP / O = a rQ / O + a rP / Q
(2.2)
Aplicando a equação de Coriolis (2.1) em (2.2),
v P / a = v Q / a + v P / b + ωb / a × rP / Q
(2.3)
A equação (2.3) é importante para os cálculos nos sistemas de navegação inercial
das aeronaves, quando se tem diferentes sistemas de coordenadas envolvidos, por
exemplo, o sistema de coordenadas do corpo da aeronave onde o sistema de
navegação está imerso, o do planeta Terra e um ponto de referência no espaço
tridimensional para onde a aeronave deve se dirigir.
Figura 2. 2: Velocidade e aceleração com frames em movimento.
A derivada da equação (2.3) fornece a aceleração de P em relação a Fa .
Analisando a equação (2.3) da esquerda para a direita, os dois primeiros termos são
definidos no frame Fa e fornecem a aceleração em relação ao mesmo. O terceiro
termo está definido em outro frame e necessita que se aplique a equação (2.1) de
Coriolis. No quarto termo, a velocidade angular torna-se vetor de aceleração angular,
tal que:
8
a P / a = aQ / a + ( a P / b + ωb / a × v P / b ) + α b / a × rP / Q + ωb / a × ( v P / b + ωb / a × rP / Q )
Reagrupando os termos, tem-se:
a P / a = a P / b + aQ / a + α b / a × rP / Q + ωb / a × (ωb / a × rP / Q ) + 2ωb / a × v P / b
aceleração
aceleração
aceleração centrípeta
relativa
de Coriolis
(2.4)
aceleração
total
aceleração de transporte de P em Fa
α : vetor de aceleração angular;
a: vetor de aceleração de translação.
As equações (2.3) e (2.4) são importantes para compor os sistemas de
navegação inerciais e a partir delas é que são possíveis estimações de variações
angulares e de translado da aeronave.
Quatérnions e vetores
A teoria dos quatérnions foi desenvolvida por W. R. Hamilton (1805-1865) na
tentativa de generalizar números complexos em espaços tridimensionais. Hoje, ela é
vastamente aplicada em simulações, robótica, navegação, controle de atitude e
computação gráfica (Stevens e Lewis; 2003) devido às vantagens que sua
representação fornece, tal como contornar problemas de singularidade matricial. Os
quatérnions são apresentados de forma bastante resumida, pois o foco deste
trabalho utiliza as representações por ângulos de Euler.
Hamilton introduziu a seguinte forma aos quatérnions:
x0 + x1i + x2 j + x3 k
e com a generalização dos complexos no espaço tridimensional
i 2 = j 2 = k 2 = ijk = −1,
ij = k , jk = i, ki = j = −ik
Quatérnions obedecem às leis da álgebra, exceto a multiplicação que é não
comutativa:
r = ( p0 + p1i + p2 j + p3k ) × ( q0 + q1i + q2 j + q3 k )
r = p0 q0 + p0 q1i + p0 q2 j + p0 q3 k + p1q0i + p1q1i 2 + ...
podendo ser escrito na forma:
⎡ r0 ⎤ ⎡ p0
⎢r ⎥ ⎢ p
⎢ 1⎥ = ⎢ 1
⎢ r2 ⎥ ⎢ p2
⎢ ⎥ ⎢
⎣ r3 ⎦ ⎣ p3
− p1
p0
p3
− p2
− p2
− p3
p0
p1
− p3 ⎤ ⎡ q0 ⎤
p2 ⎥⎥ ⎢⎢ q1 ⎥⎥
− p1 ⎥ ⎢ q2 ⎥
⎥⎢ ⎥
p0 ⎦ ⎣ q3 ⎦
9
Interpretando i, j , k como vetores unitários, os quatérnions podem ser
tratados como ( q0 + q ) , onde q é a parte vetorial do quatérnion com os componentes
q1 , q2 , q3 ao longo de i, j , k :
⎡p ⎤
p = ⎢ r0 ⎥
⎣p ⎦
⎡q ⎤
q = ⎢ 0r ⎥
⎣q ⎦
onde as componentes são tomadas em um sistema de referência r.
A multiplicação de dois quatérnions é indicada pelo operador “ ∗ ”:
p0 q0 − p ⋅ q
⎡
⎤
p∗q = ⎢
r⎥
⎣⎢( p0q + q0p + p × q ) ⎦⎥
Pelo tratamento vetorial aos quatérnions, as operações com os mesmos se
dão da mesma forma que com os vetores. Para um estudo mais completo,
aconselha-se o livro de Stevens e Lewis (2003).
Análise matricial da cinemática: transformações lineares
Considerando a equação matricial v = Au onde v e u são matrizes (n x 1) e
A , uma matriz (n x n) constante, cada elemento de v é combinação linear dos
elementos de u . Logo, v = Au é uma transformação linear da matriz u .
Supondo que em outra análise se utilize um novo grupo de variáveis através
de uma transformação linear reversível com uma matriz de transformação L , nãosingular, tal que exista uma L−1 , as novas variáveis correspondentes a u e v são:
u1 = Lu
e
v1 = Lv
A nova relação de transformação entre as variáveis é:
v1 = LAu = LAL−1u1
Matrizes
de
transformação
são
utilizadas
para
efetuar
rotação
de
coordenadas tal como as rotações com ângulos de Euler e com quatérnions, assim
como a conversão entre representações com ângulos de Euler para quatérnions,
mostrada em (2.5) e vice-versa.
q0 = ±(cos φ/2 cos θ/2 cosψ/2 + sen φ/2 sen θ/2 senψ/2)
q1 = ±(sen φ/2 cos θ/2 cosψ/2 − cos φ/2 sen θ/2 senψ/2)
q2 = ±(cos φ/2 sen θ/2 cosψ/2 + sen φ/2 cos θ/2 senψ/2)
(2.5)
q3 = ± (cos φ/2 cos θ/2 senψ/2 − sen φ/2 sen θ/2 cosψ/2)
Em (2.5) tanto faz utilizar + ou - , desde que seja único de q0 a q3 .
10
Rotações de Euler
Para melhor visualizar a orientação de aeronaves em um sistema de
coordenadas cartesianas em relação a outro, a representação de Euler é muito útil.
Descrita por três sucessivas rotações, no ramo aeronáutico as rotações de Euler são
feitas em ordem específica em torno de cada um dos três eixos cartesianos ou
planos de rotação.
A forma de efetuar estas rotações é utilizar uma matriz de transformação
chamada matriz cosseno de direção (2.6), que por questão de simplificação é
apresentada diretamente neste trabalho apenas ressaltando que os elementos
unitários e nulos correspondem à coordenadas que não mudam durante a rotação
do sistema de coordenas “a” para o “b”:
⎡ xb ⎤ ⎡ cos μ
⎢ ⎥ ⎢
⎢ yb ⎥ = ⎢ − senμ
⎢⎣ zb ⎥⎦ ⎢⎣ 0
senμ
0 ⎤ ⎡ xa ⎤
cos μ 0 ⎥⎥ ⎢⎢ ya ⎥⎥
0
1 ⎥⎦ ⎢⎣ za ⎥⎦
(2.6)
Utilizando a equação (2.6) já se pode efetuar uma rotação de coordenadas
como uma sequência de rotações de planos. De praxe, na aviação, se executa as
rotações, pela regra da mão direita e em ordem, através dos eixos z, y, x, conforme:
- Rotação no eixo z: ângulo ψ positivo (yaw ou guinada);
- Rotação no eixo y: ângulo θ positivo (pitch ou arfagem);
- Rotação no eixo x: ângulo φ positivo (roll ou rolagem).
Para melhor ilustrar, dá-se um exemplo utilizando um sistema de referência
na Terra (sistema r) e outro fixo ao “corpo” da aeronave (sistema b, do inglês, bodyfixed). Utilizando a convenção NED (do inglês, North East Down), onde o eixo z é
positivo apontado para baixo, a matriz de transformação pelas três rotações
sucessivas é obtida por:
0
0 ⎤ ⎡cos θ 0 − senθ ⎤ ⎡ cosψ senψ 0 ⎤
⎡1
⎢
u = ⎢0 cos φ senφ ⎥⎥ ⎢⎢ 0
1
0 ⎥⎥ ⎢⎢ − senψ cosψ 0 ⎥⎥ u r
⎢⎣0 − senφ cos φ ⎥⎦ ⎢⎣ senθ 0 cos θ ⎥⎦ ⎢⎣ 0
0
1⎥
⎦
b
Cb / r : matriz de transformação
Cb / r
⎡
cos θ cosψ
⎢
= ⎢( − cos φ senψ + senφ senθ cosψ )
⎢⎣ ( senφ senψ + cos φ senθ cosψ )
cos θ senψ
( cos φ cosψ + senφ senθ senψ )
( −senφ cosψ + cos φ senθ senψ )
− senθ ⎤
⎥
senφ cos θ ⎥
cos φ cos θ ⎥⎦
(2.7)
11
A matriz de transformação (2.7) possibilita encontrar outras relações com
rotações, como as que sofrem variações no tempo, podendo-se obter os
componentes de velocidade e aceleração das aeronaves. Por exemplo:
ω
b
b/r
ωbb/ r
⎡φ ⎤
⎛ ⎡0⎤
⎢ ⎥
⎜
= ⎢ 0 ⎥ + Cφ ⎜ ⎢⎢θ ⎥⎥ + Cθ
⎜ ⎢0⎥
⎢0⎥
⎝⎣ ⎦
⎣ ⎦
0
⎡ P ⎤ ⎡1
⎢
⎥
⎢
= ⎢Q ⎥ = ⎢0 cos φ
⎢⎣ R ⎥⎦ ⎢⎣0 − senφ
⎡0⎤⎞
⎢0⎥⎟
⎢ ⎥⎟
⎢⎣ψ ⎥⎦ ⎟⎠
− senθ ⎤ ⎡ φ ⎤
⎢ ⎥
senφ cos θ ⎥⎥ ⎢ θ ⎥
cos φ cos θ ⎥⎦ ⎢⎣ψ ⎥⎦
(2.8)
onde, em (2.8), P, Q e R são símbolos padrões para denotar componentes de taxas
de variação de roll, pitch e yaw respectivamente.
2.1.2 Dinâmica de corpos rígidos
Nesta seção serão destacadas algumas equações de estados importantes
para o entendimento da simulação e análise de movimento de aeronaves tal como
as que são usadas pela ferramenta de simulação Aerosim. Para um estudo mais
amplo, recomenda-se a leitura dos livros de Jan Roskam (1979), Uy-Loi Ly (1997) e
Stevens e Lewis (2003).
Movimento Angular
Usando o CM (centro de massa) da aeronave como ponto de referência, a
dinâmica rotacional pode ser separada da translacional (Stevens e Lewis; 2003).
Sendo assim, para o estudo da dinâmica rotacional, dá-se um exemplo utilizando
dois frames, Fi como referencial inercial e Fb fixo ao corpo rígido do veículo; v CM / i a
velocidade do CM no frame inercial; ωb / i a velocidade angular de Fb com respeito a
Fi e, M A,T a soma dos momentos aerodinâmicos e de propulsão sobre o CM. Seja o
vetor de momento angular no frame inercial Fi , calculado no CM de um corpo rígido
e denotado por h , sua derivada tirada no frame inercial é:
i
h = M A,T
12
Para determinar o vetor de momento angular, considere um elemento de massa δ m
com vetor posição r relativo ao CM. Sua velocidade inercial é dada por:
v = v CM / i + ωb / i × r
O momento angular desse elemento de massa em relação ao CM é igual ao
momento de translação em relação ao CM:
δ h = r × vδ m = r × v CM / iδ m + r × (ωb / i × r ) δ m
(2.9)
Integrando (2.9) sobre todos os elementos de massa, fornece:
h = ωb / i ∫ ( r ⋅ r ) dm − ∫ r ( r ⋅ ωb / i ) dm
ω
b
b/r
⎡P⎤
= ⎢⎢Q ⎥⎥
⎢⎣ R ⎥⎦
e
(2.10)
⎡x⎤
r = ⎢⎢ y ⎥⎥
⎢⎣ z ⎥⎦
b
(2.11)
Substituindo (2.11) em (2.10):
⎡ P ( y 2 + z 2 ) dm − Q xy dm − R xz dm ⎤
∫
∫
⎢ ∫
⎥
2
2
b
⎢
h = Q ∫ ( x + z ) dm − R ∫ yz dm − P ∫ yx dm ⎥
⎢
⎥
⎢ R ( x 2 + y 2 ) dm − P zx dm − Q zy dm ⎥
∫
∫
⎣ ∫
⎦
(2.12)
As várias integrais nos componentes de (2.12) são definidas como momentos
e produtos cruzados de inércia:
(
)
Momento de inércia sobre o eixo x: J xx = ∫ y 2 + z 2 dm
Produto cruzado de inércia: J xy ≡ J yx = ∫ xy dm
Substituindo estas definições em (2.12), fornece:
⎡ J xx − J xy − J xz ⎤ ⎡ P ⎤
⎢
⎥
h = ⎢ − J xy J yy − J yz ⎥ ⎢⎢Q ⎥⎥ ≡ J bωbb/ i
⎢ − J xz − J yz J zz ⎥ ⎢⎣ R ⎥⎦
⎣
⎦
b
(2.13)
matriz de inércia
A matriz J presente em (2.13) é conhecida como matriz de inércia de corpo
rígido. Ela é uma matriz real, simétrica e normalmente, por simplicidade, considerada
constante para corpos com distribuição de massa fixa, podendo ser calculada ou
determinada experimentalmente.
De posse de hb e diferenciando em relação a Fb ,
M A,T = i h = b h + ωb / i × h
(2.14)
13
M bA,T = J b bω bb/ i + Ωbb/ i J bωbb/ i
(2.15)
onde Ωbb/ i é uma matriz de produto cruzado de ωbb/ i . Rearranjando a equação (2.15)
se obtêm a equação de estados para a velocidade angular:
ω bb/ i = ( J b ) ⎡⎣ M bA,T − Ωbb/ i J bωbb/ i ⎤⎦
−1
b
(2.16)
A equação de estados (2.16) é vastamente utilizada em simulação e análise
de movimento de corpos rígidos podendo ser resolvida numericamente para ωbb/ i
quando se tem a matriz de inércia e o vetor de torque previamente calculado.
Equações de movimento de Euler
A equação de estados de velocidade angular (2.16) toma uma forma mais
simples quando a inversa da matriz de inércia for diagonal. A forma mais simples é
conhecida como equações de movimento de Euler.
Seja o vetor de torque com os componentes
M bA,T
⎡A⎤
= ⎢⎢ m ⎥⎥
⎢⎣ n ⎥⎦
as equações de Euler são:
(
P =
J y − J z ) QR
Jx
+
A
Jx
( J − J x ) RP + m
Q = z
Jy
Jy
(
R =
J x − J y ) PQ
Jz
+
(2.17)
n
Jz
As equações em (2.17) envolvem permutação cíclica dos componentes de
rotação angular e inércia, o que acarreta em um acoplamento inerente, pois,
qualquer variação angular sobre dois eixos gera uma aceleração sobre o terceiro
eixo. Este fenômeno é chamado de acoplamento inercial e será mais bem
compreendido a partir do Capítulo 3, quando algumas análises são efetuadas para
avaliar os acoplamentos do UAV Aerosonde.
14
Movimento de translação do CM
O movimento de translação de um veículo sobre a Terra pode assumir duas
formas de análise, uma mais completa, onde a rotação do planeta, achatamento dos
pólos e variações na forma de sua superfície são consideradas e outra mais simples,
que desconsidera estes fatores.
De acordo com Stevens e Lewis (2003), a análise mais completa é necessária
para simulações com precisão em veículos mais rápidos que 2000 ft/s (ou 2195
km/h) ou quando se deseja simular um sistema de navegação em longas distâncias.
Neste trabalho, o UAV utilizado nas simulações não voará a mais de 100 km/h,
descartando a necessidade da análise completa.
Para simulação de aeronaves em baixa velocidade sobrevoando uma
pequena região do planeta Terra sem a necessidade de alta precisão de posição, é
aceitável negligenciar termos centrípetos e de Coriolis das equações, assumindo
uma superfície plana e inercial. Deste modo, seja um frame Fv no veículo aéreo com
velocidade angular nula relativa a um frame de referência inercial Fe , com ωb / v ≡ ωb / e
e com o sistema geográfico de coordenadas em Fv , alinhado com um sistema plano
tangente T nas vizinhanças do veículo. A atitude da aeronave pode ser descrita por
ângulos de roll, pitch, yaw e uma matriz cosseno de direção Cb / n (corpo do veículo
com respeito ao sistema geográfico, também chamado de sistema de navegação
local que tem sua origem no CM do veículo):
Cb / n = fn ( Φ )
e
n
b
p CM
/ T = Cn / b v CM / e
= H (Φ )ωb
Φ
b/e
b
b
v CM
/e =
(2.18)
1 b
FA,T + Cb / n g n − Ωbb/ e v bCM / e
m
ω bb/ e = ( J b ) ⎡⎣M bA,T − Ωbb/ e J bωbb/ e ⎤⎦
b
−1
onde:
pCM / T : posição do CM do veículo em relação a origem do sistema T ;
v CM / e = e p CM / T : vetor de velocidade do CM em Fe ;
ωb / e : velocidade angular de Fb relativo a Fe ;
Φ : ângulos de Euler do sistema fixo ao corpo em relação ao sistema NED.
15
Levando em consideração o ar ao redor da aeronave, tem-se o vetor de
velocidade relativa:
v brel = v bCM / e − Cb / n vWn / e
(2.19)
com vW / e igual ao vetor velocidade do vento relativo a Fe .
Nesta análise simplificada, o vetor aceleração gravitacional possui apenas um
componente não nulo, g D = 9,80665 m/s2 :
gt = [0 0 g D ]
Organizando as equações de (2.18) em um vetor de estados, tem-se:
n
X T = ⎡( pCM
/T )
⎣⎢
T
ΦT
(v
) (ω )
T
b
CM / e
T
b
b/e
⎤
⎦⎥
Para se obter a equação do vetor de aceleração de translação do CM, derivase a equação (2.19) em relação ao frame fixo ao corpo da aeronave Fb eliminando
v CM / e :
b
v rel =
1 b
FA,T + g − ωb / e × v rel − e v W / e
m
(2.20)
O último termo de (2.20) pode ser utilizado para injetar perturbações no modelo, tal
como rajadas de ventos. Se igualado a zero, considerando ventos estáveis, e
introduzindo componentes do frame Fb , a equação do vetor de aceleração será:
b
v brel =
1 b
FA,T + Cb / n g n − Ωbb/ e v brel
m
(2.21)
Nesta seção verificou-se que o comportamento de um veículo aéreo rígido é
determinado pelas equações de força, momento e cinemática, tendo algumas outras
dependências mais fracas não mencionadas neste trabalho e que estão acopladas
as equações apresentadas. Algumas destas dependências não apareceram devido a
simplificações no modelo escolhido para o planeta Terra, ficando mais evidente com
o estudo da geodésia e da navegação terrestre.
2.2 Geodésia
Geodésia é uma área da matemática que estuda a forma e a área da Terra,
que é comumente descrita como um esferóide oblato, ou seja, com os pólos
levemente achatados. Essa forma se desenvolveu quando a Terra ainda estava em
16
formação a partir de uma nuvem de gases que, devido ao efeito de rotação,
concentrou maior força centrífuga na região equatorial do que na região dos pólos,
causando o achatamento ou compressão de aproximadamente 0,3% em relação ao
diâmetro equatorial (Jeppesen General Navigation Manual; 2005).
Pesquisas mais recentes, baseadas em informações de satélites, apontaram
que, na verdade, a Terra teria o formato levemente parecido ao de uma pêra, com
seu maior diâmetro, de apenas algumas dezenas de metros a mais, ao sul da linha
do equador. A questão é que para simular o movimento de um veículo aeroespacial
ao redor da Terra é necessário ter um modelo matemático que a descreva.
A partir das observações e medidas pelo estudo da geodésia, uma série de
agências definiu diferentes modelos para a Terra. Por exemplo, alguns países
europeus definiram o modelo ED50 (European Datum 1950), a França utilizou o NTF
(Nouvelle Triangulation de France) do ano de 1970 e os Estados Unidos definiu o
WGS-84 (World Geodetic System 1984), utilizado para a efeméride dos satélites do
sistema GPS, atualmente mais utilizado no mundo e que é normalmente implantado
nos diversos sistemas de navegação.
Neste trabalho, o modelo utilizado é o WGS-84, que vem integrado aos
diagramas de simulação do Aerosim no MATLAB. Como mencionado no final da
seção 2.1.2, um simples esferóide de superfície plana seria suficiente para as
simulações e projetos dos controladores, já o WGS-84 atende a simulações para
grandes velocidades sobre grandes áreas da superfície do planeta, o que é mais que
suficiente para cumprir os objetivos deste trabalho.
Além da forma, é necessário o conhecimento de um modelo de gravitação e
de rotação do planeta. A superfície equipotencial do campo gravitacional terrestre,
que coincide com o nível médio do mar (m.s.l., do inglês, mean sea level) e que se
estende continuamente abaixo dos continentes, é chamada de geoid. A distribuição
irregular de massa da Terra faz com que o geoid seja uma superfície ondulada,
conforme a Figura 2.3 apresenta.
17
Figura 2. 3: Geoid e definições de altura.
Na Figura 2.3, vê-se que a noção de vertical local é definida a partir da
direção com que a Terra atrai, por exemplo, um peso de prova, sendo normal à
superfície do geoid e o ângulo que este forma com a normal do esferóide, chama-se
deflexão vertical. O esferóide do WGS-84 tem seu centro localizado no CM da Terra
e é dividido de um em um grau (1º) tanto para a latitude como para a longitude.
Alguns parâmetros importantes deste modelo são apresentados na Figura
2.4:
Figura 2. 4: Modelo esferoidal oblato da Terra.
a ≡ 6378137, 0 m
a −b
f =
≡ 1/ 298, 257223563
a
b = 6356752 m
(a
e=
2
− b2
)
1/ 2
≈ 0, 08181919
a
ωE ≡ 7, 2921150 ×10−5 rad/s
GM ≡ 3986004, 418 × 108 m3 /s2
18
onde a, f , b, e, ωE e GM são respectivamente: eixo semi-maior, achatamento, eixo
semi-menor, excentricidade, taxa de rotação e a constante gravitacional da Terra.
2.3 Navegação terrestre
Na seção 2.1, tratou-se de frames de referência e sistemas de coordenadas
sem levar em conta as definições destes sistemas. Nesta seção, abordam-se as
nomenclaturas mais comuns dos frames e sistemas de coordenadas assim como os
sistemas de navegação e o método escolhido para efetuar a navegação do
Aerosonde durante as simulações.
2.3.1 Frames de referência e sistemas de coordenadas
Convenções NED e ENU
NED (do inglês, north, east and down): convenção para sistemas de
coordenadas que possuem seus eixos alinhados com o norte, leste e para baixo,
onde “baixo” refere-se ao vetor normal ao esferóide (ver Figura 2.3);
ENU (do inglês, east, north and up): similar ao NED mas com eixos alinhados
com o leste, norte e para cima.
Sistemas e frames
ECI (do inglês, Earth-centered inertial): sistema de coordenadas com sua
origem no CM da Terra, eixos no plano equatorial e ao longo do eixo de rotação. Um
frame de referência no ECI é denotado por Fi , sendo um frame inercial não rotativo,
mas que translada com o CM da Terra;
ECEF (do inglês, Earth-centered, Earth-fixed): sistema de coordenadas com
sua origem na superfície da Terra (sistema de plano tangente ou geográfico), eixos
no plano equatorial e ao longo do eixo de rotação. Frames no ECEF são denotados
por Fe ;
Plano-tangente: sistema geográfico com sua origem na superfície do planeta
e denotado pelo frame Ft ;
19
Vehicle-carried: sistema geográfico carregado dentro do veículo com sua
origem no CM do mesmo e aqui denotado por frames Fv ;
Vehicle body-fixed: sistema fixado ao corpo do veículo no seu CM e com
eixos alinhados conforme a convenção NED ou ENU. Frame denotado por Fb ;
Vehicle stability-axes e wind-axes: sistemas de coordenadas baseados nos
eixos de estabilidade e do vento respectivamente, sendo que o primeiro tem seus
eixos alinhados em relação ao veículo quando este encontra-se em vôo estável e
nivelado, ou seja, variações angulares nulas ( P , Q , R = 0 ). O segundo refere-se aos
ventos. Frames nestes sistemas são denotados por Fs e Fw , respectivamente.
2.3.2 Sistema sexagesimal e posição na Terra
O sistema sexagesimal utiliza o fato de que uma rotação em sentido horário a
partir do norte, leste, sul e oeste e de volta ao norte, é um círculo de 360º (Jeppesen
General Navigation Manual; 2005). O norte é definido como 000º ou 360º, o leste
como 090º, o sul 180º e o oeste 270º. Logo, as direções na Terra podem ser
medidas em graus no sentido horário a partir do norte.
Navegação é um processo fundamental quando se quer levar uma aeronave
de um ponto a outro, sendo necessário o uso de sistemas de referência de posição
sobre a superfície da Terra.
Em uma superfície plana, uma posição pode ser precisamente definida
utilizando um sistema de coordenadas cartesianas pelo uso de dois eixos
perpendiculares, x e y. Mas para superfícies esféricas, é necessário utilizar
coordenadas baseadas em ângulos, conhecidas como longitude e latitude. Definidas
de forma semelhante ao sistema cartesiano, envolvem círculos mutuamente
perpendiculares ao longo da superfície do planeta. Alguns destes círculos são
importantes destacar:
Equador: Grande Círculo cujo plano está a 90º do eixo de rotação da Terra
dividindo-a em dois hemisférios (norte e sul). Ao longo do arco do Equador estão as
direções leste e oeste;
20
Meridianos: são Semi-Grandes Círculos que unem os pólos norte e sul e
cruzam o Equador a 90º. O meridiano de Greenwich é um dos mais conhecidos,
pois é o meridiano principal que define o eixo longitudinal, análogo ao eixo y.
Figura 2. 5: Longitude e latitude em relação ao Equador e o meridiano de Greenwich.
Na Figura 2.5 ilustra-se um exemplo de como obter a posição de um ponto A
a x quilômetros ao leste do Meridiano de Greenwich e y quilômetros ao norte do
Equador ou, em termos de longitude e latitude, β º de Greenwich e α º do Equador.
Em resumo, na navegação com veículos aéreos, os conceitos de longitude e
latitude são a forma mais comum de levar a aeronave de um ponto a outro de forma
precisa, sendo para isto, necessário conhecer a posição inicial e final e ter
instrumentos que possibilitem conhecer direção, altitude, velocidade, direção dos
ventos, etc., para poder determinar a direção até o ponto desejado e o deslocamento
necessário.
2.3.3 Método de navegação adotado
Neste trabalho, um simulador de GPS é utilizado para determinar a posição
do UAV Aerosonde e assim permitir a navegação através dos waypoints, que são
pontos específicos de longitude, latitude e altitude pelos quais a aeronave deve
passar durante o vôo. O método adotado é baseado no de Randy C. Hoover (2004),
que utiliza um gerador de waypoints para que o usuário, em uma fase de pré-vôo,
21
insira uma lista de waypoints que posteriormente são utilizados pelo piloto
automático, o qual obtém a trajetória para o waypoint e monitora se o mesmo é
alcançado.
O algoritmo usado para a navegação do Aerosonde lida com o erro entre a
posição desejada e a posição atual, gerando uma trajetória longitudinal e lateral a
ser percorrida. A trajetória longitudinal é obtida simplesmente pela subtração da
altitude desejada pela atual. Já a trajetória lateral, é um pouco mais complexa e
utiliza um método por triangulação.
O método da triangulação consiste em subtrair o valor de dois vetores e
calcular a tangente inversa do resultado (Hoover; 2004), fornecendo diretamente a
direção ou heading correto para o waypoint. Os cálculos do sistema de navegação
são explicados em detalhes no Capítulo 5 juntamente com o projeto do controlador
fuzzy de Mamdani de direção.
2.4 Aerodinâmica básica
Os
modelos
matemáticos
utilizados
neste
trabalho
contêm
dados
aerodinâmicos da aeronave como um todo. No entanto, é importante compreender
como estes dados são obtidos, por partes, examinando os efeitos das forças e
momentos que agem sobre as superfícies do avião, começando pelo estudo dos
aerofólios.
2.4.1 Aerodinâmica da seção de um aerofólio
A forma de um aerofólio determina suas propriedades aerodinâmicas, e para
ilustrar os parâmetros que as afetam, a Figura 2.6 apresenta uma seção de um
aerofólio em duas dimensões e supostamente de comprimento infinito.
22
Figura 2. 6: Seção de um aerofólio e principais parâmetros.
A Figura 2.6 ilustra um aerofólio submetido a uma corrente de ar ao redor de
sua superfície com linhas tangentes a ela garantindo um fluxo estável e constante. A
linha de corda (chord line) é a linha de referência para descrever a forma do
aerofólio, definindo-o como simétrico ou assimétrico. A linha média ou linha de
abaulamento (camber line) divide a superfície superior e inferior em distâncias
iguais, definindo a espessura do aerofólio. Estes parâmetros combinados definem a
faixa de velocidade da aeronave.
Imaginando um eixo perpendicular ao plano da Figura 2.6, que passa através
da linha de corda e inclinando o aerofólio pela rotação nesse eixo, forma-se o ângulo
α , conhecido como ângulo de ataque. Com a variação deste ângulo, duas forças
aerodinâmicas perpendiculares são alteradas, o empuxo (que garante a sustentação
e elevação) e o arrasto devido à fricção do ar com a superfície do aerofólio,
conforme visto na Figura 2.6. O ângulo de ataque muitas vezes é confundido com o
ângulo de pitch θ , no entanto, a diferença é que o ângulo de pitch é formado pela
rotação no eixo y em relação ao sistema de referência da aeronave, do tipo bodyfixed, como mencionado na seção 2.3.1, ou seja, é o ângulo entre o eixo longitudinal
da aeronave e o horizonte. Já o ângulo de ataque, é a inclinação em relação ao fluxo
de ar ou à direção do vento (Wind-axes). Um exemplo clássico para visualizar essa
diferença é o caso de um avião em aproximação para pouso. Normalmente, o avião
se aproxima com o “nariz” levantado, por exemplo, com pitch θ = 10º , mas está
descendo, o que é um fator aerodinâmico, logo, dependente do ângulo de ataque α .
Este, só coincide com θ quando a força resultante no eixo vertical z for nula.
Na análise lateral também ocorre certa confusão entre o ângulo de guinada
(yaw) ψ e o ângulo de derrapagem (sideslip) β , sendo estes respectivamente
análogos aos casos de θ e α . A Figura 2.7 apresenta tais diferenças.
23
Na Figura 2.7, além de se definir os ângulos aerodinâmicos, também são
apresentados os eixos relativos ao vento incidente e o eixo de estabilidade, sendo
esse último, o que define a direção do vôo (ou flight path). Na mesma figura, as três
principais superfícies móveis de controle de um avião (sendo estas as únicas
utilizadas neste trabalho) são ilustradas, e aqui apresentadas suas funções:
Figura 2. 7: Definição das superfícies de controle, de ângulos e eixos aerodinâmicos.
Ailerons: superfícies móveis localizadas nas asas que atuam simetricamente
e são responsáveis pela rotação no eixo longitudinal x. Efetuam o movimento de
rolagem (roll) denotado pelo ângulo φ ou ângulo de bank (bank angle);
Profundor (elevator): superfície móvel localizada na cauda da aeronave no
estabilizador horizontal responsável pela rotação no eixo lateral y. Efetua o
movimento de arfagem (pitch) denotado pelo ângulo θ ;
Leme (rudder): superfície móvel localizada na cauda da aeronave no
estabilizador vertical responsável pela rotação no eixo z. Efetua o movimento de
guinada (yaw) denotado pelo ângulo ψ de direção ou heading.
24
Além destas superfícies móveis, pode-se citar os flaps, slats, trens de pouso,
etc, que por restrições de tempo, não são utilizadas neste trabalho e, para mais
informações, recomenda-se a leitura da referência de Jan Roskam (1985). Outro
comando comum em aviões motorizados é o acelerador (ou throttle), que varia a
potência entregue ao motor da aeronave.
A convenção de sinais para as superfícies de controle adotadas neste
trabalho utiliza a convenção comum da indústria aeronáutica (Figura 2.7 e Tabela
2.1):
Tabela 2. 1: Convenções de sinais para as superfícies de controle aerodinâmicas.
Deflexão
Sinal
Efeito
Elevator
Para baixo
Positivo
Baixa o nariz da aeronave
Rudder
Para a esquerda
Positivo
Guinada do nariz para a esquerda
Ailerons
Aileron direito baixo Positivo
Asa esquerda desce e a direita sobe
Com o que foi abordado até o presente momento já se tem informação
suficiente para seguir ao estudo da planta de controle e entender quais variáveis
estão envolvidas nas simulações. Os principais ângulos e vetores conhecidos neste
capítulo irão compor o vetor de estados que será derivado a partir do modelo não
linear do Aerosonde, que é o veículo aéreo não tripulado apresentado no Capítulo 3.
25
Capítulo 3
Aerosim 1.2, introdução aos UAVs e ao Aerosonde
Neste capítulo, abordam-se as funcionalidades do Aerosim, sigla para
Aeronautical Simulation Blockset (ou, em português, Conjunto de blocos para
Simulações Aeronáuticas), e suas principais características, faz-se uma introdução
aos UAVs destacando os tipos e suas funções, finalizando-se com a apresentação
da planta de controle, o Aerosonde.
3.1 Aerosim 1.2
O Aerosim encontra-se na versão 1.2 sendo um conjunto de ferramentas
desenvolvidas para simulação de modelos dinâmicos não-lineares de aeronaves
com 6-DOF. Desenvolvido e distribuído pela empresa Unmanned Dynamics
(www.u-dynamics.com), que tem por objetivo prover softwares para UAVs e foi
fundada em 2002 pelo engenheiro de controle de aeronaves, Marius Niculescu. O
Aerosim é distribuído gratuitamente para fins educacionais no site da empresa.
Essa ferramenta é composta de diversos blocos para o Simulink/MATLAB, da
empresa MathWorks, e engloba uma biblioteca de modelos completos de aviões
personalizáveis através da modificação de parâmetros. Dentre esses modelos, o
modelo não-linear do Aerosonde está presente juntamente com alguns exemplos de
simulações para facilitar o aprendizado e rapidamente familiarizar o usuário com a
ferramenta.
Além dos modelos de aeronaves, o Aerosim possibilita a imersão em um
ambiente de simulação de uma atmosfera com ou sem ventos constantes, bem
como simulação de rajadas de ventos e turbulências, e utiliza o modelo da Terra
baseado no WGS-84. Proporciona meios de efetuar transformações para vários
frames de referência, como body-axes, wind-axes, frames geográficos, ECI e ECEF;
e blocos de conversão de unidades comuns utilizadas em engenharia aeronáutica.
Podendo simular toda a dinâmica da aeronave assim como o ambiente ao
redor dela, o Aerosim também conta com uma maneira de reduzir a abstração da
simulação indo além das análises gráficas comuns, podendo exibir animações
26
tridimensionais ao conectar o Simulink/MATLAB ao simulador de vôo Microsoft Flight
Simulator, distribuído pela empresa Microsoft, que fica responsável somente em
exibir a simulação que ocorre no MATLAB. Pode-se, também, conectar um
controlador de jogo tipo joystick para pilotar manualmente a aeronave através do
Aerosim.
O Aerosim requer a versão do MATLAB 6, ou mais recente. Sua instalação e
ligação
do
mesmo
ao
simulador
Microsoft
Flight
Simulator
é
explicada
detalhadamente no APÊNDICE A.
3.1.1 Interface
O Aerosim inclui ao Simulink as bibliotecas de blocos mostradas na Figura
3.1. Atenção especial é dada ao item de aviões completos (Complete Aircraft) onde
se encontra o principal bloco utilizado neste trabalho, o de modelo de avião com 6DOF.
Figura 3. 1: Interface da biblioteca de blocos do Aerosim 1.2 no Simulink.
Na Figura 3.1, o bloco responsável pela simulação do modelo não-linear de
avião com 6-DOF possui três entradas e quinze saídas para possibilitar uma
infinidade de testes. Este bloco contém sub-blocos do pacote Aerosim interligando
27
os diversos modelos: do avião, da atmosfera e do geoid da Terra. Os canais de
entrada/saída mais comuns e utilizados neste trabalho são:
Controls: vetor de entrada 7x1 de controles do avião como: flaps, profundor,
ailerons, leme, acelerador, mistura e ignição. Os controles aerodinâmicos estão em
radianos, o acelerador vai de 0 a 1, a mistura é dada pela relação
ar / fluxo combustível e ignição 0 ou 1 para desligado e ligado respectivamente;
Winds: vetor de entrada 3x1 das componentes de velocidade dos ventos
dadas em um frame de referência geográfico do tipo NED;
RST: entrada booleana para efetuar reset nos estados da aeronave;
States: saída com o vetor 15x1 de estados do avião, sendo as variáveis em
ordem: posição (latitude, longitude e altitude), velocidade (u, v, w, em m/s), atitude
(q0, q1, q2, q3, do quatérnion), variações angulares ( P, Q, R , em rad/s), massa de
combustível (kg) e rotação do motor (RPM). A atitude em quatérnions pode ser
diretamente convertida para ângulos de Euler utilizando os blocos de transformação
que acompanham o Aerosim;
Sensors: vetor 18x1 das saídas dos sensores, como: posição com GPS,
velocidade
com
GPS,
acelerômetros,
giroscópios,
dados
atmosféricos
e
magnetômetro;
VelW: vetor 3x1 de saída da velocidade do avião no sistema de referência
wind-axes, contendo: velocidade do fluxo de ar, ângulo de derrapagem β e ângulo
de ataque α ;
Euler: saída do vetor 3x1 dos ângulos de Euler: roll, pitch e yaw;
MSL: altitude do avião em relação ao nível do mar;
AGL: altitude do avião acima do nível do solo (AGL, do inglês, Above Ground
Level).
O bloco da Figura 3.1 é uma generalização de um modelo com 6-DOF, sendo
necessário alimentá-lo com parâmetros da aeronave e suas condições iniciais.
Dando um duplo-clique com o mouse sobre o bloco, a tela de parâmetros se abre no
Simulink. É importante ressaltar a função de alguns itens:
Aircraft configuration file: arquivo de configuração da aeronave. Deve-se
fornecer o caminho para o arquivo do tipo MAT do MATLAB que contém os dados da
aeronave que são pertinentes a montagem e obtenção do seu modelo não linear. No
caso deste trabalho, será utilizado o arquivo aerosondecfg.mat, pré-configurado com
28
os dados do UAV Aerosonde e mostrado no ANEXO A. Para construção de arquivos
MAT personalizados, o Aerosim trás um script em formato m-file do MATLAB que
“pergunta” ao usuário alguns dados da aeronave para gerar o arquivo de
configuração. Para maiores informações, deve-se consultar o manual do usuário do
Aerosim (www.u-dynamics.com);
Initial position: vetor 3x1 de posição inicial da aeronave dado em latitude,
longitude e altitude;
Initial Velocities: vetor 3x1 de velocidade inicial da aeronave dado em
componentes NED;
Initial attitude: vetor 4x1 de atitude inicial dada em quatérnions. Pode-se
utilizar o bloco de transformações do Aerosim para converter atitude em ângulos de
Euler para quatérnions;
Initial angular rates: vetor 3x1 de variações angulares iniciais ( P, Q, R );
Initial fuel mass: massa inicial da quantidade de combustível à bordo da
aeronave (dada em quilogramas);
Initial engine speed: velocidade inicial de rotação do motor (em RPM);
Sample time: período de amostragem (em segundos).
Com as informações apresentadas nesta seção pode-se ter uma idéia das
principais funções disponíveis no Aerosim para promover as simulações deste
trabalho. Tendo a ferramenta sido apresentada, passa-se ao estudo das aplicações
em que simulações com veículos aéreos são utilizadas.
3.2 Introdução aos UAVs
Os UAVs podem ser divididos em três categorias, os que voam de forma
autônoma, os que são pilotados remotamente (RPVs, do inglês, Remotely Piloted
Vehicle) e os que dependem de pilotagem remota em alguns trechos da missão, por
exemplo, na decolagem e no pouso, mas seguem autônomos nos demais trechos.
Normalmente, um sistema de UAV é composto de um veículo aéreo, uma ou
mais estações de controle e monitoramento em solo, carga útil (doravante designada
como payload) e um link de dados para troca de informações ar/solo (Introduction to
UAV System; 2002).
29
Nenhuma operação com UAVs teria propósito se não fossem pelos payloads.
Normalmente são os subsistemas mais caros dentro de um sistema de um UAV.
Geralmente incorporado por câmeras para alta e baixa luminosidade, sistema de
reconhecimento de alvo, marcadores a laser, radares, dispositivos de sensoriamento
de variáveis ambientais e químicas, armas, etc.
UAVs ganharam reconhecimento global pela sua integração efetiva em
operações militares de defesa e ataque. Devido a seus motores de baixo ruído e
difícil detecção visual ou por radares, tornaram-se elemento essencial nos campos
de batalha. Além destas aplicações, os UAVs têm potencial em operações não
militares,
como:
“no
campo
de
monitoração
ambiental,
climatológica
e
biodiversidade; mapeamento; prospecção topográfica, mineral e arqueológica;
levantamento de ocupação urbana e em áreas rurais; inspeção de linhas de
transmissão, oleodutos e gasodutos; observação de obras e represas; monitoração
de estradas e rios; geração e retransmissão de sinais de comunicação; operação em
zonas de perigo ou de desastre; vigilância, busca e salvamento” (Ramos e Bueno;
2007).
Estes veículos podem variar quanto à estrutura, tamanho e modelo, sendo os
tipos mais comuns os aviões, helicópteros e dirigíveis. Alguns países podem ser
citados como pioneiros na pesquisa e desenvolvimento de UAVs, como: Estados
Unidos, Israel, Austrália e países europeus. No âmbito nacional, desde o início da
década de 1980, começou-se a desenvolver pesquisas nessa área, conduzidas no
Centro Técnico Aeroespacial (CTA). Desde então, diversos projetos passaram a ser
desenvolvidos até os dias atuais, como exemplos: projeto Helix (1992-1995) para um
helicóptero não tripulado da empresa Gyron Tecnologia com parceria de diversas
universidades brasileiras e, culminando em 1996, no desenvolvimento do projeto
AURORA (Autonomous Unmanned Remote Monitoring Robotic Airship), para um
dirigível não tripulado (Ramos e Bueno; 2007); projeto ARARA (1999), sigla para
Aeronaves de Reconhecimento Assistidas por Rádio e Autônomas, iniciado pela
USP de São Carlos, utilizando aviões (Neris; 2001); projeto SiDeVAAN (2004), sigla
para Simulação e Desenvolvimento de Veículos Aéreos Autônomos e Nãotripulados, realizado pela Universidade Federal de Minas Gerais (UFMG) com aviões
(Campos et al.; 2007).
30
3.3 Aerosonde UAV
O Aerosonde UAV é a aeronave escolhida para ser utilizada como planta a
ser controlada neste trabalho. Isto se deve ao fato de que a ferramenta Aerosim já
trás um modelo não-linear em blocos do Simulink para esta aeronave. O Aerosonde
é desenvolvido e operado globalmente há mais de dez anos pelas empresas
Aerosonde
Pty
Ltd
(AePL)
e
Aerosonde
North
America
(AeNA)
(www.aerosonde.com).
Este UAV foi desenvolvido com o intuito de operar em vôos econômicos,
longos e a partir de qualquer região do planeta. Possibilitando vôos de até 32 horas,
foi o primeiro UAV a cruzar o Oceano Atlântico em 1998 e a executar vôos dentro de
furacões em 2005 e 2007 para estudo e predição dos mesmos. Com
aproximadamente 3 metros de envergadura de asa e 13 quilogramas, o Aerosonde
utiliza sistema de navegação por GPS e uma estação de monitoramento e comando
em solo para ser operado. Permite payloads de até 2,3 quilogramas, velocidades de
cruzeiro que variam de 60 km/h a 110 km/h, razão de subida de até 9 km/h e
altitudes de até 6000 metros (20000 ft).
Na Figura 3.2, são apresentadas duas imagens do Aerosonde, sendo a da
esquerda uma fotografia real e a da direita, o modelo gráfico tridimensional da
aeronave no simulador Microsoft Flight Simulator.
Figura 3. 2: Fotografia e modelo gráfico tridimensional do Aerosonde UAV.
31
3.3.1 Estudo da dinâmica do Aerosonde
Para o estudo da dinâmica do Aerosonde, uma análise em malha aberta
fornece informações básicas para avaliar a estabilidade do sistema e como o mesmo
se comporta quando submetido a pequenas perturbações. Na Figura 3.3 apresentase o diagrama de blocos de simulação no Simulink, do modelo não-linear do
Aerosonde. O mesmo é considerado como um modelo tipo caixa preta apenas para
testar seu comportamento em malha aberta e obter algumas conclusões sobre as
respostas das simulações no contexto de estabilidade, para introduzir o conceito de
“trimagem” (trim) que possibilitará a obtenção de modelos linearizados em pontos de
operação específicos.
Figura 3. 3: Diagrama de simulação no Simulink do Aerosonde tipo caixa preta.
Executando a simulação do modelo da Figura 3.3, com a aeronave em malha
aberta, durante 1 minuto, sem ventos e monitorando a dinâmica lateral do
Aerosonde, ou seja, coletando as respostas para o ângulo de bank φ e de direção
(heading) ψ , quando nas entradas (Controls) o vetor [-0.1 0 0 0.4] = [profundor
ailerons leme acelerador] é aplicado, as respostas da Figura 3.4 e 3.5 são obtidas,
mostrando a instabilidade na dinâmica lateral da planta mesmo com as superfícies
de controle, ailerons e leme, estando niveladas e fixadas na posição neutra.
32
Figura 3. 4: Resposta temporal do ângulo de bank.
Figura 3. 5: Resposta temporal de heading.
A instabilidade na dinâmica lateral observada nas figuras 3.4 e 3.5 deve-se ao
torque do motor que gera um movimento espiral em torno do eixo longitudinal x. Em
33
detalhes, o ângulo de bank aumenta demasiadamente e o avião entra em “parafuso”
dando diversas voltas completas de 0º a 359º de direção.
Por tentativa e erro, pode-se chegar a um valor de entrada que estabilize as
asas do avião mantendo-as niveladas pelo menos por um curto período de tempo. O
mesmo pode ser encontrado para a dinâmica longitudinal, procurando um ângulo de
pitch que possibilite ao avião voar reto e nivelado, ou seja, com a resultante de todas
as forças que atuam sobre a aeronave, nulas. A busca por tais valores é chamada
de ajuste do trim, que dada uma velocidade, altitude e peso da aeronave, “sintoniza”
as superfícies de controle e potência do motor para que o vôo seja reto e nivelado.
O ajuste do trim (ou “trimagem”) é fundamental para que seja possível obter
modelos lineares do Aerosonde em pontos de operação específicos de velocidade,
altitude, massa e assim derivar um modelo dentro de uma região de estabilidade.
Por tentativa e erro, o valor de -0,0073 rad de deflexão dos ailerons manteve o
ângulo de bank mais próximo de zero, conforme a Figura 3.6. No entanto, por
tentativa e erro levaria muito tempo até conseguir um arranjo de todas as superfícies
de controle e aceleração que garantisse um vôo reto e nivelado. Como solução,
métodos iterativos de otimização, por exemplo, utilizando o algoritmo Simplex
(Stevens e Lewis; 2003), ou mesmo a pilotagem manual, conseguem obter
resultados satisfatórios.
Figura 3. 6: Resposta temporal de
φ
para entrada nos ailerons de -0,0073 rad.
34
Na versão 1.2 do Aerosim, um programa para o ajuste do trim e obtenção dos
modelos lineares vem acompanhado de um diagrama de simulação para o Simulink,
onde, através da função linmod do Matlab, extrai-se um modelo linear no espaço de
estados. Este programa é apresentado no ANEXO B deste trabalho. Em resumo,
este programa solicita ao usuário a velocidade (m/s), a altitude (m), o ângulo de bank
(graus) e a massa de combustível (kg).
O programa de “trimagem” solicita a informação de ângulo de bank pois podese querer analisar a dinâmica da aeronave enquanto efetua curvas estáveis de
ângulos constantes, na literatura específica comumente chamadas de steady-turns.
Quando a dinâmica do avião é analisada nas condições de vôo reto e
nivelado ou em steady-turn, pode-se dividir o estudo da dinâmica em longitudinal e
lateral, pois o acoplamento destas dinâmicas é mínimo e pode ser desconsiderado
(Stevens e Lewis; 2003). Este tipo de análise facilita na visualização e entendimento
do problema assim como no projeto dos controladores que podem ser feitos
separadamente.
Para estudar os modos ou pólos mais importantes na dinâmica do Aerosonde,
um modelo linear nas seguintes condições é obtido: velocidade = 23 m/s, altitude =
1000 m, φ =0º e 2 kg de combustível. O resultado do ajuste para este ponto de
operação é:
Entradas:
Saídas:
Elevator = -0,1429 rad
Airspeed = 23,00 m/s
Aileron = -0,0086 rad
Sideslip = 0,02º
Rudder = -0,0010 rad
AOA = 4,32º
Throttle = 0,5698 (~57%)
Bank = -0,02º
Pitch = 4,32º
Heading = 0,04º
Altitude = 1000,00 m
Estados:
u = 22,93 m/s (componente no eixo x)
v = 0,01 m/s (componente no eixo y)
w = 1,73 m/s (componente no eixo z)
p = -0,00 º/s, q = -0,00 º/s, r = -0,00 º/s
35
φ = -0,02º, θ = 4,32º, ψ = 0,04º
Motor = 4835 RPM
Utilizando os valores do vetor de entrada obtidos pelo programa de
“trimagem” na simulação da planta não linear, pode-se perceber variações menores
nos ângulos φ , θ e na velocidade, que se mantém em torno de 23 m/s conforme as
figuras 3.7, 3.8 e 3.9.
Figura 3. 7: Resposta temporal de
φ
do modelo não-linear com os valores de entrada trimados.
Se estas simulações fossem feitas com o modelo linearizado, utilizando os
mesmos valores de entrada, mesmas condições iniciais e sem perturbações, a
aeronave permaneceria no estado de equilíbrio.
Apesar de se ter utilizado os valores obtidos pela “trimagem”, ainda assim, a
planta não-linear sofre os efeitos do torque do motor com o avanço do tempo, no
entanto, de forma menos agressiva do que sem os ajustes. De qualquer modo,
percebe-se que o Aerosonde não voa reto e nivelado sem um controlador humano
ou automático.
36
Figura 3. 8: Resposta temporal da velocidade relativa (airspeed) para o ponto de operação em 23 m/s.
Figura 3. 9: Resposta temporal do ângulo
θ
para o ponto de operação de trimagem.
Os modelos lineares longitudinal e lateral obtidos nesta trimagem são
apresentados e estudados na seção seguinte.
37
3.3.1.1 Dinâmica longitudinal
O modelo longitudinal linear obtido a partir da trimagem para um vôo do
Aerosonde a 23 m/s, 1000 m de altitude, asas niveladas e 2 kg de massa de
combustível foi:
⎡ -0, 2029 0, 6110 -1, 7044 -9, 7991 -0, 0001 0, 0100 ⎤
⎢ -0,5670 -3,8086 22, 4291 -0, 7397 0, 0010
0 ⎥⎥
⎢
⎢ 0, 4906 -4, 2213 -4,3901
0
0, 0000 -0, 0078 ⎥
A=⎢
⎥
0
1, 0000
0
0
0 ⎥
⎢ 0
⎢ 0, 0753 -0,9972
0
22,9997
0
0 ⎥
⎢
⎥
0
0
-0, 0379 -2, 6886 ⎦
⎣ 29, 7630 2, 2466
0
⎡ 0,3132
⎤
⎢ -1,9847
⎥
0
⎢
⎥
⎢-27,5486
⎥
0
B=⎢
⎥
0
0
⎢
⎥
⎢
⎥
0
0
⎢
⎥
0
349, 2659 ⎦
⎣
;
⎡ 0,9972 0, 0753 0 0 0
⎢-0, 0033 0, 0434 0 0 0
⎢
C=⎢ 0
0
1 0 0
⎢
0
0 1 0
⎢ 0
⎢⎣ 0
0
0 0 1
0⎤
0 ⎥⎥
0⎥
⎥
0⎥
0 ⎥⎦
onde,
Vetor de estados: x = [u w q theta h Omega]
Vetor de entrada: u = [elevator throttle]
Vetor de saída: y = [Va alpha q theta h]
u = componente da velocidade (m/s) no eixo x;
w = componente da velocidade (m/s) no eixo z;
q = variação angular (rad/s) de θ ;
h = altitude (m);
Omega = rotação do motor (RPM);
Va = airstream velocity ou airspeed (velocidade relativa m/s).
Com as asas niveladas, o acoplamento da dinâmica longitudinal e lateral pode
ser desprezado. Desta forma pode-se estudar os modos naturais envolvidos em
cada dinâmica separadamente, fornecendo ao projetista maior conhecimento sobre
a planta, como é mostrado nos exemplos a seguir:
38
Simulação em malha aberta aplicando pulsos de elevator/profundor
Modo natural: short-period
Utilizando condições iniciais nulas é aplicado na entrada da planta um pulso
duplo no elevator, de amplitude 0,1 rad de 1 a 1,5 segundos e de -0,1 rad de 1,5 a 2
segundos. Este pulso duplo é chamado de doublet e possui média nula na intenção
de restaurar as condições de vôo quando os pulsos finalizam. A Figura 3.10 mostra
a resposta temporal do ângulo de pitch e do ângulo de ataque após o doublet no
elevator. As respostas das duas saídas não coincidem plenamente em forma nem
em duração ao distúrbio que as causou. Estas respostas são características da
aeronave e representam o modo natural denominado short-period, no qual α e θ
variam juntos (ver seção 2.4.1), causando pouca mudança no ângulo γ = θ − α ,
conhecido como flight path angle ou, em português, ângulo de translado do vôo.
Na Figura 3.11, inspecionam-se outras variáveis envolvidas, como a altitude e
velocidade, que pela ordem de grandeza em que são normalmente utilizadas, podem
ser consideradas praticamente constantes, sendo então o short-period responsável
por variações consideráveis apenas em α e θ .
Figura 3. 10: Resposta temporal de
α eθ
ao doublet no elevator e modo natural short-period estimulado.
39
Figura 3. 11: Resposta temporal de h e Va ao doublet no elevator.
Modo natural: phugoid
Fazendo outra análise sobre as respostas das figuras 3.10 e 3.11 observa-se
outro modo natural sendo estimulado, o modo phugoid. Quando os efeitos do shortperiod começam a cessar, ainda sobra uma pequena amplitude de oscilação pouco
amortecida que representa a troca de energia potencial e cinética enquanto a
aeronave varia sua atitude. Este efeito é melhor visualizado na Figura 3.11, onde
percebe-se que quando a velocidade aumenta a altitude diminui, o que indica que o
ângulo de translado γ se tornou negativo e a aeronave começou a descer e ganhar
velocidade. Devido ao aumento do fluxo de ar passando pelos aerofólios da
aeronave, ele aumenta sua sustentação invertendo o movimento e ganhando altitude
e perdendo velocidade. Este ciclo se repete até as oscilações sumirem e a aeronave
retornar ao seu estado inicial.
Simulação em malha aberta aplicando-se doublet na entrada de aceleração
Aplicando-se a mesma técnica do exemplo anterior, só que desta vez, na
segunda entrada, referente ao acelerador, com pulsos de 0,5 de 1 a 4 segundos e
de -0,5 de 4 a 7 segundos, tem-se a resposta para α e θ na Figura 3.12:
40
Figura 3. 12: Resposta temporal de
α eθ
ao doublet no throttle e oscilação em
θ
devido ao phugoid.
Examinando a Figura 3.12, percebe-se que o ângulo de ataque é pouco
afetado enquanto que o de pitch sofre grandes variações. Conclui-se então, que o
short-period é pouco afetado, no entanto, as oscilações do phugoid acompanham a
atitude de pitch durante toda a simulação e, analisando a velocidade e altitude na
Figura 3.13, observa-se que a altitude e velocidade variam junto com θ . Logo, a
entrada do acelerador influencia predominantemente no modo phugoid.
Análise dos pólos e zeros de transmissão do sistema
Pela análise dos pólos do sistema e com o auxílio das figuras 3.10 a 3.13 é
possível identificar os pólos do sistema que estão envolvidos com os modos shortperiod e phugoid. Este estudo é importante para verificar se a planta possui pólos
instáveis e também para determinar a freqüência e o amortecimento. Isto dará mais
informações ao projetista sobre as dificuldades que poderão ser encontradas quando
for desenvolver os controladores. Dependendo do fator de amortecimento e da
freqüência desses pólos, pode-se avaliar a qualidade do vôo, baseado tanto em
pilotagem manual como em automática, sendo alguns casos tão críticos que a
pilotagem manual só poderia ser feita se houvesse a presença de um sistema de
41
controle de estabilização entre o piloto humano e a planta. Este tipo de análise e
sistema será abordado no Capítulo 4.
Figura 3. 13: Resposta temporal de h e Va ao doublet no throttle. Phugoid estimulado.
Os pólos, amortecimentos, freqüências e constantes de tempo do modelo
longitudinal são mostrados na Tabela 3.1:
Tabela 3. 1: Pólos, fatores de amortecimento, freqüências e constantes de tempo do modelo longitudinal.
Pólo
Amortecimento
−6, 49 ×10−4
1
−5,35 ×10−2 ± j 5, 73 ×10−1
9,30 ×10−2
−2, 79
1
−4,10 ± j 9, 77
0,387
Freqüência (rad/s)
Cte. Tempo (s)
1540
5, 76 ×10−1
phugoid
short − period
0,3584
10, 6
Os pólos envolvidos com a dinâmica do short-period e phugoid são sempre
complexos conjugados, conforme destacados na Tabela 3.1. Isto se deve ao que já
fora mostrado nas figuras 3.10 a 3.13, onde a dinâmica de ambos os modos era
42
oscilatória. A grande diferença está na freqüência e no amortecimento. O shortperiod geralmente tem a maior freqüência (menor período) e maior amortecimento.
No modo phugoid dá-se o contrário. Estes resultados não são animadores, pois se
percebe que apesar de estável, o Aerosonde possui um modo short-period subamortecido e de período muito curto. O modo phugoid também se mostra
problemático, já que as oscilações realmente levam muito tempo para cessar devido
ao baixíssimo fator de amortecimento. Em compensação, tem um período maior que
facilita a pilotagem. O sistema também possui dois zeros de transmissão localizados
no semiplano esquerdo do plano-s, sendo então o sistema de fase mínima com
z1 = -3, 4766 e z2 = 0, 0003 . Perceba que z2 está muito próximo do primeiro pólo da
Tabela 3.1, podendo-se fazer um cancelamento destes.
As informações obtidas no estudo da dinâmica longitudinal já possibilitaram
prever algumas dificuldades que deverão ser levadas em consideração quando os
controladores forem projetados, tal como problemas de amortecimento e forte
acoplamento entre as variáveis. Para completar o estudo da dinâmica do Aerosonde,
dá-se seqüência a análise da dinâmica lateral.
3.3.1.2 Dinâmica lateral / direcional
O modelo linear da dinâmica lateral obtido foi:
⎡-0,5895 1, 7312 -22,9345 9, 7991
⎢ -3,8720 -19, 0490 9,1681
0
⎢
A = ⎢ 0, 6278 -2, 4709 -0,9582
0
⎢
1, 0000
0, 0755
0
⎢ 0
⎢⎣ 0
0
1, 0028
0
2,9486 ⎤
⎡ -1,1552
⎢-101, 4284 1,8250 ⎥
⎢
⎥
B = ⎢ -3,9992 -18, 6309 ⎥
⎢
⎥
0
0
⎢
⎥
⎢⎣
⎥⎦
0
0
;
⎡0, 0435
⎢ 0
⎢
C=⎢ 0
⎢
⎢ 0
⎢⎣ 0
0⎤
0 ⎥⎥
0⎥
⎥
0⎥
0 ⎥⎦
0
1
0
0
0
0
0
1
0
0
0
0
0
1
0
0⎤
0 ⎥⎥
0⎥
⎥
0⎥
1 ⎥⎦
onde,
Vetor de estados: x = [v p r phi psi]
43
Vetor de entradas: u = [aileron rudder]
Vetor de saídas: y = [beta p r phi psi]
v = componente da velocidade (m/s) no eixo y;
p = variação angular (rad/s) de φ ;
r = variação angular (rad/s) de ψ .
Da mesma maneira como foi feito o estudo da dinâmica longitudinal, serão
feitas algumas simulações em malha aberta (MA) na forma de exemplos com o
modelo lateral.
Simulação MA com aplicação de doublet no rudder
A aplicação de um doublet na segunda entrada, a do rudder, excita os três
modos naturais que se deseja observar no comportamento dinâmico lateral do
modelo. Os pulsos aplicados têm amplitude 0,5 rad de 1 a 2 segundos e de -0,5 rad
de 2 a 3 segundos de simulação. As figuras 3.14 a 3.16 apresentam as respostas
relacionando as variáveis β , φ , ψ , de sideslip, roll e heading, respectivamente.
Figura 3. 14: Resposta de
β
e
φ
ao doublet aplicado no rudder.
44
Figura 3. 15: Resposta de
β
e ψ ao doublet aplicado no rudder.
Figura 3. 16: Resposta de
φ eψ
ao doublet aplicado no rudder.
45
Nas figuras 3.14 e 3.15 observa-se que quando o rudder é estimulado
ocorrem variações no ângulo de guinada ψ e também no ângulo de rolagem φ . Este
efeito é decorrente do modo natural de oscilação denominado dutch roll. Logo, o
rudder da aeronave produz movimento tanto de guinada como de rolagem. Percebe-
se também, que ocorre derrapagem devido à variação no ângulo de sideslip β ,
sempre contrária ao movimento de φ e ψ . Para ilustrar melhor este efeito, toma-se
como exemplo um veículo terrestre. Um carro, ao iniciar uma curva, faz o passageiro
ser "puxado" para o lado contrário ao da curva devido à atuação de forças
centrífugas. Normalmente o carro não sofre a derrapagem, pois o atrito com o asfalto
é muito grande, mas com o avião, esta derrapagem é percebida e considerável
conforme as figuras 3.14 e 3.15. A derrapagem cessa após alguns segundos quando
o estímulo dos pulsos também cessa.
A Figura 3.16 compara as variações que ocorrem nos ângulos φ e ψ ,
servindo para avaliar outro modo natural denominado roll subsidence. Este modo é
responsável pela taxa de variação do ângulo φ e é normalmente estável. Fornece
uma idéia da máxima taxa de rolagem da aeronave.
O terceiro modo natural é designado por spiral, que tem efeito semelhante ao
estudado no início do capítulo, provocado pelo torque do motor e que causou o
movimento espiral fazendo a aeronave entrar em “parafuso”. Apesar desta
comparação, este modo normalmente conta com grandes constantes de tempo, o
que muitas vezes a própria estrutura da aeronave é projetada para eliminar os
efeitos espirais ou o piloto, automático ou manual, facilmente o eliminam. Para
facilitar o entendimento, dá-se um exemplo de um carro desalinhado, que necessita
de constantes correções na trajetória, por parte do piloto, para que o veículo se
mova em linha reta.
No Aerosonde, o caso da própria estrutura eliminar o modo espiral não
ocorre, como pode ser visto, na Figura 3.16, φ e ψ progressivamente desviando-se
da posição neutra após o pulso ter cessado.
Análise dos pólos e zeros de transmissão do sistema
O modelo linear lateral do Aerosonde possui dois zeros de transmissão,
ambos no semiplano esquerdo do plano-s, sendo o sistema de fase mínima com
z1 = -148, 0479 e z2 na origem. Quanto aos pólos do sistema, um deles encontra-se
46
na origem e se cancela com z2 . Os demais pólos, mostrados na Tabela 3.2, são
relativos aos três modos naturais envolvidos na dinâmica lateral. O dutch roll é subamortecido e de curto período, o que prejudica a pilotagem principalmente na fase
de pouso sob rajadas de ventos. O modo espiral é instável e, portanto, precisa ser
modificado por algum sistema de controle ou estabilização, que é o foco dos
sistemas que serão estudados e projetados no Capítulo 4.
Tabela 3. 2: Pólos, fatores de amortecimento, freqüências e constantes de tempo da dinâmica lateral.
Pólo
Amortecimento
Freqüência (rad/s)
Cte. Tempo (s)
−1, 22 ± j 5,39
0, 221
5,53
−18, 2
1
0, 0549
6, 46 ×10−2
−1
15, 4798
dutch roll
roll subsidence
spiral
47
Capítulo 4
Visão geral dos sistemas de controle de aeronaves
Os modos naturais estudados no Capítulo 3 podem ser divididos em
diferentes categorias, por exemplo: os que envolvem principalmente movimentos
rotacionais, como o short-period, roll subsidence e dutch roll; os modos que
envolvem movimentos translacionais e de dinâmica muito mais lenta, como o
phugoid e, por último, o modo spiral, dependente de ambos os movimentos, masque
é normalmente excitado apenas por forças aerodinâmicas fracas. Todos esses
modos influenciam nas respostas da aeronave.
A boa resposta de uma aeronave aos comandos de manobra é determinada
em parte pela velocidade dos modos rotacionais. As freqüências destes modos e o
pouco amortecimento podem tornar difícil ou impossível para um piloto controlar a
aeronave. Para solucionar tais situações, tal como verificado no Aerosonde no
Capítulo 3, com a presença de pólos instáveis e pouco amortecidos, é necessário
utilizar algum sistema de controle automático que modifique os modos naturais para
torná-los viáveis à pilotagem. Estes tipos de sistemas são conhecidos como SAS (do
inglês, Stability Augmentation Systems), responsáveis por melhorar a estabilidade da
aeronave. Se estes sistemas, além de aumentarem a estabilidade, também
garantirem tipos particulares de respostas às entradas de controle de um piloto
(automático ou manual), são denominados de CAS (do inglês, Control Augmentation
Systems), podendo alterar as características dos modos naturais para diferentes
características de vôo.
Os sistemas de piloto automático são relacionados normalmente às funções
“superiores” de controle e podem utilizar CAS e SAS em camadas inferiores. Um
exemplo de piloto automático poderia ser um segurador de atitude de pitch (Pitch
attitude hold), utilizando em uma camada mais baixa na hierarquia da malha de
controle, um CAS para controlar as variações das taxas de pitch (por exemplo, Q
em rad/s) e um SAS que aumentasse o amortecimento do modo short-period.
A Tabela 4.1 apresenta tipos comuns de SAS, CAS e funções de pilotos
automáticos (Stevens e Lewis; 2003). Os tipos foram deixados em inglês para
familiarização do leitor com a linguagem normalmente utilizada tanto no Brasil como
48
no exterior, onde dampers são amortecedores, rates as taxas de variação e holders
são seguradores:
Tabela 4. 1: Tipos comuns de sistemas de aumento de estabilidade, de controle e pilotos automáticos.
SAS
CAS
Roll damper
Roll rate
Pitch damper
Pitch rate
Yaw damper
Lateral / direcional
Pilotos Automáticos
Pitch attitude hold
Altitude hold
Speed hold
Roll-angle hold
Para garantir o projeto de controladores satisfatórios, algum critério de
comparação deve ser utilizado como requisito de desempenho. Na literatura,
normalmente são encontrados alguns critérios de comparação baseados em
opiniões de operadores humanos para veículos tripulados, o que não exclui sua
utilização em UAVs, já que estes podem necessitar de intervenção humana através
da pilotagem remota. De qualquer modo, o foco deste trabalho é a pilotagem
autônoma de veículos que geralmente transportam payloads sensíveis a trepidações
e solavancos. Portanto, seguir os requisitos de qualidade e desempenho das
aeronaves tripuladas aumenta as chances de se obter bons resultados.
Diversas
pesquisas
no
ramo
aeroespacial
são
voltadas
a
essas
especificações de qualidade de vôo e vários métodos têm sido propostos. A
princípio, tais estudos parecem simples, mas, na prática, ainda não é possível
especificar, com precisão, critérios de projeto de sistemas de controle para modificar
a dinâmica de aeronaves (Stevens e Lewis; 2003). Alguns desses métodos
consistem em analisar pólos e zeros, ou resposta em frequência, e comparar com
especificações conhecidas e opiniões de pilotos.
4.1 Especificações de pólos/modos
As especificações apresentadas neste trabalho foram extraídas de Stevens e
Lewis (2003) e são baseadas nas Especificações Militares de Qualidade de Vôo dos
Estados Unidos. Elas são definidas de acordo com a classe de aeronave, fases de
vôo e por níveis de qualidade de vôo. Neste trabalho, são abordadas apenas as
49
especificações resumidas para as classes de aviões pequenos e leves, fases de vôo
não terminais (decolagem, aproximação e pouso) e qualidade nível 1, que seria a
qualidade adequada à fase de vôo.
Especificações para os modos:
Phugoid: ζ ≥ 0, 04 ;
Short-period: ζ ≥ 0,30 ;
Roll subsidence: constante de tempo máxima τ = 1, 4 s ;
Spiral: modo pode ser instável desde que não dobre em amplitude após 20 s;
Dutch Roll: ζ ≥ 0, 08 .
4.2 Sistemas de aumento de estabilidade (SAS)
A maioria das aeronaves, tripuladas ou não, requerem algum tipo de sistema
que melhore sua estabilidade. Existem aviões que são tão inerentemente instáveis
que seria praticamente impossível o vôo destes sem um sistema automático de
controle de estabilidade. Os SAS normalmente utilizam as medições de sensores
como giroscópios para a tomada das variações das taxas angulares da aeronave.
Essas medidas são processadas e realimentadas aos servomecanismos que atuam
nas superfícies de controle aerodinâmicas. Deste modo, um momento aerodinâmico
proporcional à velocidade angular e à aceleração angular pode ser gerado e produzir
um efeito de amortecimento nos movimentos da aeronave.
Convencionalmente, os SAS podem ser projetados separadamente para a
dinâmica longitudinal e para a lateral-direcional. Isto é possível devido ao
desacoplamento destas dinâmicas em boa parte das fases de vôo, tal como fora
mostrado no Capítulo 3. Nos casos em que existe grande influência entre as
dinâmicas longitudinal e lateral, como em curvas com grandes variações no ângulo
de bank, os controladores deverão garantir o desacoplamento (ver seção 4.3).
Os controladores que serão desenvolvidos para aumentar a estabilidade e
controlar o Aerosonde UAV, neste trabalho, são apontados como sub-controladores
dentro de uma hierarquia de comando, utilizando sistemas SAS e CAS em uma
única estrutura. No entanto, para facilitar o entendimento dos projetos, é importante
50
ressaltar de forma independente alguns tipos de sistemas SAS e CAS para conhecer
quais as variáveis envolvidas e o que se deseja com cada um destes sistemas. Os
sistemas SAS mais comuns são então descritos detalhadamente na forma de
exemplos.
4.2.1 Pitch SAS
A principal função de um sistema do tipo Pitch SAS é garantir amortecimento
e freqüência natural satisfatória para o modo natural short-period estudado no
Capítulo 3. Foi visto que este modo natural envolve tanto o ângulo de ataque α
como as variações no ângulo de pitch (e também no pitch rate, q ). Normalmente,
faz-se realimentação apenas do pitch rate se o modo natural short-period possuir
amortecimento adequado de acordo com as especificações de qualidade. Caso
contrário, fecha-se outra malha de realimentação com a saída do ângulo de ataque
ou do próprio ângulo de pitch, conforme mostrado na Figura 4.1:
Figura 4. 1: Estrutura de um Pitch SAS com realimentação do ângulo de ataque e pitch rate.
Onde: u é o sinal de controle; ue o erro; δ e é a deflexão do pofundor/elevator;
kα e kq são os ganhos do controlador.
O sistema apresentado na Figura 4.1 realimenta a saída do ângulo de ataque
e quando as perturbações são verificadas nessa variável a realimentação leva esta
informação ao atuador que procura gerar um ângulo de pitch restaurador.
Malhas externas de controle geralmente são fechadas ao redor do Pitch SAS
para promover funções de sistemas CAS e de Piloto Automático.
51
4.2.2 SAS Lateral-Direcional ou Yaw Damper
O SAS que aumenta a estabilidade da dinâmica lateral das aeronaves é
chamado de Yaw Damper e está presente na maioria dos sistemas de controle de
vôo de aviões comerciais e militares, na intenção de modificar os modos naturais
roll-subsidence e dutch-roll pela realimentação das variações do ângulo de bank
(roll, φ ) aos ailerons, e das variações do ângulo de guinada (yaw, ψ ) ao leme. O
Yaw Damper procura eliminar os efeitos de derrapagens durante curvas ou quando a
aeronave é submetida a perturbações como rajadas de ventos laterais. A estrutura
do Yaw Damper é apresentada na Figura 4.2:
Figura 4. 2: Estrutura de um Yaw Damper.
Onde
r1
e
r2
são as referências para as variações angulares
p
e
r
respectivamente; ua é o comando de aileron e ur o de leme/rudder; δ a é a deflexão
dos ailerons; δ r é a deflexão do leme/rudder; k p e kr os ganhos do controlador.
Na Figura 4.2 é apresentado um bloco designado como “Filtro Washout” que
tem como entrada o yaw rate ( r ) e saída rw . Este filtro não é essencial, mas
normalmente está presente nos aviões manualmente pilotados. Isto se deve ao fato
de que o Yaw Damper tenta se opor a atuação, por exemplo, de um comando do
piloto ao leme, fazendo com que o piloto tenha que aplicar comando acima do
normal para obter o efeito de guinada. Sendo assim, o Filtro Washout é um filtro
passa-alta que facilita a pilotagem manual, logo, quando a aeronave encontra-se em
uma condição de regime, com pouca ou nenhuma variação em r , o Yaw Damper
52
procura restaurar a atitude quando submetida a pequenas perturbações. Para
perturbações muito grandes, como a da entrada de um comando de leme enviado
pelo piloto, o filtro modifica as características do Yaw Damper para possibilitar a
manobra. Neste trabalho, o Filtro de Washout não será implementado, mas sim um
controlador multivariável que procura ter as características do Yaw Damper,
modificadas a partir de um índice de desempenho, como um CAS-SAS em um só.
4.3 Sistemas de aumento de controle (CAS)
Enquanto que os sistemas SAS procuram modificar as características
dinâmicas das aeronaves no sentido de dar-lhes melhores condições de pilotagem e
estabilidade, os sistemas CAS procuram controlar certas variáveis para possibilitar
que estas possam desempenhar tarefas específicas, como por exemplo, rastrear
com precisão um percurso ou um alvo de forma automática. Além desses, imagine
que em aeronaves de alto desempenho, como caças, que precisam executar
manobras rápidas e precisas e geralmente a grandes velocidades, controlar as
variáveis de variações angulares previne estresse da estrutura da aeronave e
desmaio do piloto devido aos efeitos provocados pela força g (em termos de n vezes
a aceleração da gravidade na Terra) elevada.
Exemplos comuns de sistemas CAS são os controladores de roll rate, pitch
rate e controle lateral-direcional. O exemplo com aviões tipo caças militares é um
caso extremo, mas analisando o caso do Aerosonde UAV e outras aeronaves lentas
e não acrobáticas, ainda assim é importante a presença dos sistemas CAS para
eliminar acoplamentos entre a dinâmica longitudinal e lateral em curvas
coordenadas, reduzir o efeito de forças indesejáveis e garantir melhor qualidade de
vôo.
Quando um avião inicia uma curva, seja ela pela atuação dos ailerons ou pelo
leme, conforme visto no Capítulo 3, geralmente ocorre rotação em torno do eixo
longitudinal, ou seja, rolagem. Sabe-se também, que a sustentação da aeronave dáse pela atuação de forças no sentido de baixo para cima, sendo boa parte desta
sustentação proveniente das forças que atuam nas asas. Quando um avião executa
um roll, mesmo que de pequeno ângulo de bank, é como se o comprimento das asas
houvesse encurtado, reduzindo a sustentação conforme mostrado na Figura 4.3:
53
Figura 4. 3: Perda de sustentação durante o roll influenciando na dinâmica longitudinal e lateral.
Além da perda de sustentação, percebe-se que no sistema de referência body-axis o
que era ângulo de pitch passa a ter certa influência no ângulo de guinada e viceversa. Logo, para reduzir tal acoplamento, deve-se controlar as variações angulares
para eliminar as forças laterais (e longitudinais) indesejáveis geradas pela rolagem e
evitar mudança na atitude da aeronave. Se cada controlador executar corretamente
sua função, os efeitos negativos da manobra serão eliminados. Para visualizar
melhor os sistemas CAS, serão apresentadas algumas estruturas de controle que
facilitarão o entendimento dos mesmos antes destes serem implementados com as
técnicas de controle robusto multivariável.
4.3.1 Pitch-Rate CAS
Na Figura 4.4 é apresentado um diagrama de blocos representativo de um
possível sistema CAS. Veja que este se assemelha ao Pitch SAS da seção 4.2.1
exceto pelo controlador C ( s ) presente na entrada do SAS para regular o valor de q
em uma referência desejada. Normalmente, como citado na seção anterior, se quer
regular q em zero e com erro de regime nulo para que, por exemplo, durante uma
curva coordenada, não haja alteração no ângulo de pitch.
54
Figura 4. 4: Diagrama de um possível pitch-rate CAS.
O controlador C ( s) poderia ser, por exemplo, um compensador Proporcional
Integral (PI) para garantir erro de regime nulo.
Estruturas de controle semelhantes poderiam ser construídas para controlar
as variações angulares dos ângulos de bank e guinada, não sendo necessário seguir
exatamente nesta mesma linha em que simplesmente se fechou uma malha externa
ao redor do sistema SAS. Diferentes estruturas podem ser experimentadas e
recomenda-se a leitura das referências de Stevens e Lewis (2003) e Roskam (1979)
para um estudo mais detalhado das estruturas de controle comumente utilizadas na
indústria aeroespacial.
4.4 Pilotos Automáticos
Diferentemente do que foi abordado sobre os sistemas SAS e CAS, que são
voltados às características de resposta dinâmica da aeronave, os pilotos
automáticos devem ser projetados para atender as especificações de erro de regime
e rejeição de distúrbios, exercendo a função de aliviar o trabalho de um piloto
humano (dentro da aeronave ou remoto) em manter uma determinada condição de
vôo por longos períodos.
Existem diversos tipos de pilotos automáticos normalmente utilizados em
aviões, como: de velocidade, altitude, atitude, etc. Alguns são tão complexos e
precisos que englobam uma infinidade de subsistemas todos interligados e
supervisionados para serem ativados e desativados no momento certo. No entanto,
na sua forma mais simples, um piloto automático de um avião geralmente utiliza os
sistemas apresentados na Tabela 4.1, incluindo os sistemas SAS e CAS em suas
malhas internas de controle.
55
4.4.1 Pitch-Attitude Hold
O Pitch-Attitude Hold é normalmente utilizado quando a aeronave está com as
asas niveladas (ângulo de bank φ = 0 ). A variável controlada é a do ângulo de pitch,
θ = γ + α , onde γ é o ângulo de direção de translação do vôo no eixo longitudinal.
Este tipo de piloto automático é utilizado para manter o nariz do avião em um ângulo
desejado e pode ser utilizado como um sub-controlador para um piloto automático de
altitude. Apenas ressaltando o que já fora comentado no Capítulo 2, somente com a
variação de θ , não significa que a aeronave irá subir ou descer, pois quem vai
determinar isso é o ângulo γ , que também depende de α . O diagrama de blocos
que descreve este piloto automático é mostrado na Figura 4.5.
Figura 4. 5: Diagrama de blocos de um piloto automático tipo Pitch-Attitude Hold.
4.4.2 Altitude Hold
O piloto automático de altitude possibilita que a aeronave seja mantida em
uma altitude fixa desejada. Nos aviões modernos de passageiros, ele normalmente
possui uma tolerância de erro de regime de ±200 ft ( ±60 m). Nos sistemas mais
modernos, além de fixar a aeronave na altitude desejada, também a conduz de uma
altitude a outra sob condições especificadas de razão de subida e limites de ângulo
de pitch.
Na Figura 4.6 é apresentado o diagrama de blocos de um piloto automático
de altitude que utiliza o controlador de ângulo de pitch da seção anterior e os
sistemas SAS e CAS nas malhas internas, possibilitando o projeto do piloto
automático de altitude com limitador de ângulo de pitch.
56
Figura 4. 6: Diagrama de blocos de um piloto automático de altitude com sub-controlador de ângulo de pitch.
4.4.3 Speed Hold
O piloto automático de controle de velocidade utiliza a medida da velocidade
do fluxo de ar que passa através da aeronave que é realimentada e comparada à
velocidade desejada.
Devido ao forte acoplamento entre a velocidade e a variável do ângulo de
pitch, o controlador de velocidade, quando em união com o controlador de altitude,
torna-se bastante complexo, pois a estrutura do Speed Hold comanda tanto o
acelerador como o profundor da aeronave para regular a velocidade. Isto influencia
diretamente no controle da altitude. Uma possível estrutura do Speed Hold seria o
fechamento de uma malha externa no Pitch-Attitude-Hold para a variação de θ no
intuito de ganhar ou perder velocidade. Um exemplo seria a estrutura do diagrama
da Figura 4.7 onde o controlador de altitude também foi incluído.
Figura 4. 7: Diagrama de blocos de um Speed Hold e Altitude Hold na mesma malha de controle.
A estrutura de controle da Figura 4.7 utiliza as duas variáveis de entrada de
controle longitudinal da planta, o acelerador/throttle e o profundor/elevator. Esta
configuração é mais complexa, pois tanto o controlador de altitude quanto o
57
controlador de velocidade afetam o ângulo de pitch. No entanto, as técnicas de
controle robusto multivariável que serão utilizadas na elaboração do Speed Hold do
Aerosonde UAV lidam com este problema de forma simples e sem a preocupação
com a estrutura do controlador.
4.4.4 Roll Angle Hold
A forma mais simples de um piloto automático de Roll Angle Hold é um
nivelador de asas que regula o ângulo de bank φ em zero pela realimentação dessa
variável aos comandos de entrada dos ailerons. O caso mais completo possibilita o
rastreamento de diferentes ângulos de bank desejados e a estrutura é mais
complexa, pois exige CAS e SAS presentes tanto na dinâmica lateral como na
longitudinal para controlar o pitch-rate ( Q ) e a derrapagem ( β ) e possibilitar curvas
coordenadas (steady-turns), tal como mencionado na seção 4.3.
O Roll Angle Hold geralmente é utilizado como sub-controlador para outros
sistemas, como por exemplo, quando está na malha interna do piloto automático de
direção que garante o vôo em um ângulo específico da bússola, tal como pode ser
visto no projeto do controlador de direção do Aerosonde UAV no Capítulo 5.
A estrutura comumente utilizada neste piloto automático é semelhante a da
Figura 4.5, apenas substituindo θ por φ e pitch-rate ( Q ou q ) por roll-rate ( P ou
p ). Com esta última estrutura é finalizado o estudo sobre os sistemas de controle de
vôo normalmente presentes nas aeronaves e que possibilitará uma melhor
compreensão do que será abordado no projeto dos controladores deste trabalho.
4.5 Técnicas de controle robusto multivariável
As técnicas de controle moderno são baseadas diretamente no modelo por
variáveis de estados tal como o que fora obtido no Capítulo 3, e ao contrário das
técnicas de controle clássico e seus sucessivos fechamentos das malhas de controle
com obtenção individual dos ganhos dos controladores, no controle moderno os
problemas são resolvidos pela solução de equações matriciais que permitem que os
58
ganhos sejam computados simultaneamente, fechando todas as malhas de controle
de uma só vez. Desta forma o projetista conta com maior liberdade e simplificação
quando o sistema possui mais de um par entrada/saída ou múltiplas malhas de
controle.
Em resumo, pela utilização de técnicas de controle moderno multivariável, os
projetos por tentativa e erro e o fechamento sucessivo das malhas de controle são
reduzidos, restando ao projetista apenas trabalhar com a seleção de um critério de
desempenho que será apresentado mais adiante neste capítulo.
Os problemas que serão abordados nesta seção envolvem regulação e
rastreamento de uma referência utilizando o método LQG/LTR. O problema de
regulação baseia-se em levar uma dada saída do sistema à zero garantindo uma
característica de resposta temporal desejada, tal como o que é exigido pelos
sistemas SAS, CAS e de pilotos automáticos. O problema de rastreamento exige
que a saída “siga” uma referência desejada, podendo esta ser uma constante não
nula ou uma referência variante no tempo.
4.5.1 Robustez e estabilidade
No projeto de sistemas de controle de aeronaves baseados em modelos
matemáticos lineares é importante ressaltar, que tais modelos, são aproximações do
modelo não linear dentro de uma região de estabilidade ou num certo ponto de
operação, seja ele estável ou não. Logo, dinâmicas de alta freqüência podem estar
sendo negligenciadas no modelo linear. Além disso, quando o avião se afasta da
região de estabilidade, ou da condição de vôo onde fora linearizado, seu
comportamento varia. Esta variação de parâmetros é um efeito de baixa freqüência
que pode tornar o sistema instável. Para compensar estas variações e erros de
modelagem, deve-se utilizar controladores que garantam a robustez e estabilidade
do sistema levando em consideração os erros de modelagem dos modos de alta
freqüência e variações de parâmetros.
Além dos erros de modelagem, os sistemas de controle de aviões devem
contornar problemas de perturbações, como, rajadas de ventos e ruídos de medição
dos sensores, garantindo desempenho satisfatório em termos de sobre-sinal, tempo
de acomodação, etc.
59
O desempenho dos controladores em termos de robustez e estabilidade no
LQG/LTR (Doyle e Stein; 1979) é garantido pela minimização de um índice de
desempenho (PI, do inglês, performance índex), ou função custo, e pela análise do
sistema no domínio da freqüência, descrita pelo diagrama de magnitude de Bode
para sistemas MIMO (do inglês, Multiple-Input Multiple-Output). Para isto, é
necessária uma representação especial para a magnitude de funções matriciais, que
são os ganhos principais inferiores e superiores, ou valores singulares (SVs, do
inglês, singular values) do sistema. Esses devem respeitar margens que garantam
estabilidade e robustez de forma análoga a utilizada no controle clássico, sendo que
para o caso MIMO, o máximo SV deve ter baixa magnitude nas altas freqüências
para garantir a estabilidade devido aos erros de modelagem e ruídos de medição,
enquanto que o mínimo SV deve ter grande magnitude nas baixas freqüências para
garantir a robustez. Estas noções ficarão mais claras na forma de exemplos que
serão apresentados nos projetos dos controladores do Aerosonde UAV.
4.5.2 Especificações de desempenho no domínio da freqüência
Resumidamente, as especificações de desempenho nas baixas freqüências
exigem que o ganho principal inferior, ou mínimo SV ( σ ), possua grande magnitude,
enquanto que o ganho principal superior, ou máximo SV ( σ ), nas altas freqüências,
possua pequena magnitude, fornecendo assim uma largura de banda ampla, porém
restrita, com frequência de corte ( ωc ) definida para que o ganho de malha GK ( jω )
tenha valor próximo de 1, ou 0 dB, em ωc . Logo, o máximo valor singular na
frequência de corte especificada deve satisfazer a equação 4.1:
σ (GK ( jωc )) = 1
(4.1)
E para garantir a robustez e precisão estática do sistema sujeito a perturbações nas
baixas freqüências ( ωd ), a especificação em termos do mínimo valor singular para
GK ( jω ) é:
σ (GK ( jω )) 1 para ω ≤ ωd
(4.2)
De praxe, para garantir que as perturbações de baixa frequência sejam atenuadas
de forma satisfatória, deve-se garantir que σ (GK ( jω )) seja maior que 40 dB.
60
Outra especificação de baixa frequência refere-se ao erro de regime nulo.
Para garanti-lo, pelo teorema do valor final (Stevens e Lewis; 2003) obtém-se a
equação 4.3, onde r é a magnitude da entrada de um degrau unitário e δ ∞ é o maior
valor que o erro de regime pode assumir.
σ (GK (0)) >
r
δ∞
(4.3)
Nas altas freqüências, a partir de certa frequência ωn , geralmente estão
presentes ruídos originados pelos sensores e dinâmicas não modeladas, logo, a
margem que garante a robustez e estabilidade é limitada pela equação 4.4:
σ (GK ( jω )) 1 para ω ≥ ωn
(4.4)
Para garantir que os ruídos de medição sejam atenuados satisfatoriamente,
de praxe, deve-se garantir que o ganho GK em (4.4) seja menor que -20 dB.
Estas margens servirão como base para o processo de loop-shaping no
projeto dos controladores LQG/LTR, onde o projetista procura definir o formato dos
SVs mínimo e máximo no diagrama de magnitude de Bode de sistemas MIMO.
4.5.3 LQG
O Regulador Linear Ótimo Quadrático (LQR) e o Filtro de Kalman podem ser
utilizados em conjunto para o projeto de um regulador dinâmico como o LQG. Uma
das características interessantes do uso do LQG é que a estrutura do compensador
é fornecida pelo processo de projeto e não necessita ser conhecida de antemão. Isto
torna o LQG uma ótima opção para o controle de sistemas muito complexos e com
um grande número de variáveis, tal como o Aerosonde UAV, que possui ordem
bastante elevada e grande número de entradas e saídas.
Para o projeto do LQR, um exemplo básico seria um sistema descrito no
espaço de estados pelas equações 4.5 e 4.6, com x(t ) ∈ R n e u (t ) o sinal de entrada
de controle por realimentação total de estados, mostrado em (4.7), onde r pode ser
uma referência a ser rastreada e K é o ganho do controlador LQR:
x = Ax + Bu
(4.5)
y = Cx
(4.6)
u = − Kx + r
(4.7)
61
Fornecendo assim o sistema em malha fechada na equação 4.8:
x = ( A − BK ) x + Br
(4.8)
O problema é que o LQR é baseado na realimentação total de estados e geralmente
nem todas as variáveis estão disponíveis para serem realimentadas. É aí que entra a
solução pelo regulador LQG, que utilizando um Filtro de Kalman, estima as variáveis
não medidas e as realimenta ao LQR. Fornecendo, assim, a equação 4.9, onde x̂ é
o vetor de estados estimados e L o ganho do Filtro de Kalman:
xˆ = ( A − LC ) xˆ + Bu + Ly
(4.9)
A nova lei de controle, substituindo os estados estimados, em (4.7), passa a ser:
u = − Kxˆ + r
(4.10)
sendo a estrutura do LQG apresentada na Figura 4.8.
A estrutura do controlador LQG faz com que os pólos de malha fechada do
sistema comportem-se identicamente como se os estados tivessem sido
realimentados sem estimação ao controlador LQR, podendo assim usufruir dos
benefícios deste regulador, como estabilidade garantida de malha fechada. No
entanto, quando o LQR e o Filtro de Kalman são utilizados nesta configuração, o
LQR “perde” algumas das suas propriedades vantajosas, tal como margem de ganho
infinita e margem de fase de 60º. Para contornar este problema, o LTR (Doyle e
Stein; 1979) é utilizado para recuperar tais propriedades do LQR e moldar por loopshaping a forma dos ganhos principais em função da freqüência e atender às
margens de robustez e estabilidade estudadas na seção 4.5.2. Para verificar
individualmente onde o LQR e o Filtro de Kalman entram em conflito, eles serão
apresentados em mais detalhes.
Figura 4. 8: Estrutura do controlador LQG.
62
4.5.3.1 LQR
Um sistema (estabilizável) em malha fechada com LQR é garantidamente
estável para o caso de realimentação total de estados. No entanto, geralmente, nem
todos os estados de um sistema podem ser medidos para satisfazer este requisito,
mas, como foi visto antes, com o uso do LQG pode-se projetar um estimador de
estados e usufruir das vantagens do regulador LQR.
Equação Algébrica de Riccati e ganho do LQR
Para se obter o ganho K ótimo do LQR de realimentação de estados em
(4.7), que minimiza o funcional em (4.11), é necessário resolver a Equação Algébrica
de Riccati (ARE, do inglês, Algebraic Riccati Equation) em (4.12) e obter a matriz
simétrica P , ( P = PT ), semidefinida positiva. De posse de P , esta é substituída na
equação 4.13 para se determinar o ganho ótimo do controlador LQR.
J=
1 ∞ T
x Qx + u T Ru ) dt
(
∫
0
2
(4.11)
0 = AT P + PA + Q − PBR −1 BT P
(4.12)
K = R −1 BT P
(4.13)
onde Q ≥ 0 e R > 0 são matrizes de pesos conhecidas que determinam o índice de
desempenho. As magnitudes relativas de Q e R devem ser selecionadas, como
parâmetros de projeto, para priorizar relativamente os requisitos de minimização dos
estados e do esforço de controle. Uma matriz R de maior magnitude fará com que o
sinal u (t ) seja minimizado, logo, otimiza-se o esforço de controle, mantendo
Ru (t )
próximo de zero. Diz-se que R muito grande penaliza mais o controle, fazendo-o
menor relativo ao vetor de estados. Por outro lado, para fazer com que x(t ) tenda a
zero mais rapidamente, seleciona-se Q com maior magnitude. Estas propriedades
fazem de Q e R os principais parâmetros de escolha no projeto do regulador LQR,
pois estes afetam diretamente na posição dos pólos de malha fechada do sistema e
nas características de resposta temporal e em frequência. Outros aspectos
importantes referem-se à parte de análise.
Na parte de análise, deve-se verificar se o sistema é detectável e
estabilizável. Se
(
)
Q , A for detectável, então ( A, B ) é estabilizável se e somente se
63
existir uma única solução semidefinida positiva de P para a equação de Ricatti em
(4.12) (Anderson e Moore; 1989). Se estes requisitos forem atendidos, o LQR
garante a estabilidade assintótica do sistema em malha fechada. Sendo a
detectabilidade uma condição menos rigorosa de observabilidade, onde se verifica
apenas se os pólos instáveis são observáveis. Logo, para garantir a detectabilidade,
todos os pólos instáveis da planta devem aparecer no índice de desempenho em
(4.11). Já a estabilizabilidade, corresponde a uma condição menos rigorosa de
controlabilidade, exigindo que apenas os pólos instáveis sejam controláveis,
garantindo-se que o controle u (t ) tenha influência suficiente sobre eles.
Na literatura específica existem diversas formas propostas para a escolha de
Q e R . No entanto, neste trabalho, tais matrizes são determinadas tomando como
base o conhecimento sobre a planta e suas variáveis de estado, utilizando matrizes
de peso modificadas, tais como as abordadas em Stevens e Lewis (2003). Isto se
deve aos acoplamentos existentes entre as variáveis e a necessidade de deixar
algumas delas “mais soltas” que outras, restringindo então a matriz do ganho.
Para um estudo mais completo sobre o assunto, recomenda-se a leitura das
referências de Doyle e Stein (1979), Anderson e Moore (1989), Green e Limebeer
(1995) e Stevens e Lewis (2003).
O foco deste estudo sobre LQR é sua aplicação, juntamente com o Filtro de
Kalman para compor a estrutura do LQG. Deve-se ressaltar que o regulador LQR
propõe-se a garantir as especificações de baixas freqüências, tal como mostrado em
(4.2) para garantir margens de ganho elevadas nas baixas freqüências. Já o Filtro de
Kalman, abordado na próxima seção, procura atenuar os efeitos das altas
freqüências, entrando em conflito com o LQR, fazendo diminuir as altas margens de
ganho do sistema.
4.5.3.2 Filtro de Kalman
Nas seções anteriores foi explicado que para poder utilizar a realimentação
total de estados com o LQR, muitas vezes, há necessidade de se estimar as
variáveis a partir das que são realmente medidas pelos sensores da planta. Além
64
disto, há necessidade de se filtrar os efeitos indesejáveis das dinâmicas de alta
frequência, erros de modelagem e ruídos de medição para garantir a estabilidade.
No ramo da aviação, o Filtro de Kalman é vastamente utilizado devido ao
elevado número de variáveis presentes, com algumas sendo muito difíceis de se
tomar medidas precisas a partir de sensores. Um exemplo é o ângulo de ataque, que
normalmente é estimado por um Filtro de Kalman a partir do pitch-rate ( q ), ângulo
de pitch ( θ ) e da aceleração normal da aeronave.
Como este trabalho é centrado somente em simulação, todas as variáveis de
estado estão disponíveis para realimentação, porém, por questões didáticas e até
mesmo pensando em uma futura aplicação em um protótipo real, a estrutura com o
Filtro de Kalman foi escolhida.
Utilizando o mesmo sistema descrito em (4.5) e (4.6), seja xˆ (t ) o estado
estimado de x(t ) , um sistema dinâmico observador de estados é descrito como:
xˆ = Axˆ + Bu + L ( y − Cxˆ )
(4.14)
xˆ = ( A − LC ) xˆ + Bu + Ly
(4.15)
ou
Pela equação 4.15 vê-se que o observador de estados é um sistema com
duas entradas, u (t ) e y (t ) , que são respectivamente ação de controle e saída da
planta, ambas conhecidas. As saídas do observador de estados são ŷ e x̂ :
yˆ = Cxˆ
(4.16)
y = y − yˆ
(4.17)
x = x − xˆ
(4.18)
Sejam y e x dados por:
onde y é o erro de estimação da saída e x o erro de estimação dos estados.
Diferenciando-se (4.18), e usando (4.5) e (4.15), tem-se que:
x = ( A − LC ) x
(4.19)
O erro de estimação inicial é dado por x ( 0 ) = x ( 0 ) − xˆ ( 0 ) , com xˆ ( 0 ) sendo o valor
inicial estimado, que geralmente é escolhido igual a zero.
A estrutura do observador de estados, tal como a utilizada com o Filtro de
Kalman, é apresentada na Figura 4.9:
65
Figura 4. 9: Estrutura da planta com observador de estados.
É necessário que o erro de estimação desapareça com o tempo para
qualquer x ( 0 ) . Para que isto ocorra,
( A − LC )
deve ser assintoticamente estável.
Logo, desde que seja obtido um ganho L que garanta esta estabilidade, os estados
serão estimados corretamente. Além do fator de estabilidade, é necessário que o
erro de estimação desapareça rapidamente, o que pode ser feito com facilidade
desde que o estimador seja assintoticamente estável e tenha dinâmica mais rápida
que a planta.
Para encontrar um ganho L tal que ( A − LC ) seja assintoticamente estável,
pode-se utilizar o problema de realimentação de estados de forma similar ao LQR
para que L aloque os pólos do sistema de acordo com as propriedades desejadas.
Então, de forma análoga ao ganho do LQR, K , em (4.13), obtido pela solução da
Equação Algébrica de Riccati em (4.12), formula-se o problema dual para L :
Substituir A, B, K
por
AT , C T , LT
LT = R −1CP → L = PC T R −1
(4.20)
0 = AP + PAT + Q − PC T R −1CP
(4.21)
A equação 4.21 é conhecida como Equação Algébrica de Riccati de Observador.
Outro aspecto desta dualidade está entre a condição de alcançabilidade1 de
( A, B )
e a observabilidade de ( BT , AT ) , com o dual para o caso do observador para
1
Um dado estado x1 é dito alcançável se existir uma entrada capaz de conduzir o sistema de um estado inicial x0
a x1 dentro de um intervalo de tempo finito.
66
a observabilidade de
( C , A)
e alcançabilidade de
( A, Q ) . Então, o sistema em
(4.19) obtido com o ganho L , e P (única e definida positiva), é assintoticamente
estável.
A grande vantagem desta dualidade é que os parâmetros de projeto
continuam sendo as definições das matrizes Q e R tal como com o LQR.
Até o presente momento, o problema foi atacado apenas do ponto de vista da
estimação dos estados. No entanto, outra vantagem do Filtro de Kalman é lidar com
as incertezas do modelo, perturbações e ruídos que atuam sobre a planta, utilizando
as técnicas de análise no domínio da frequência através dos ganhos principais.
O Filtro de Kalman baseia-se em um tratamento probabilístico do processo e
ruídos de medição, sendo caracterizado como um filtro passa-baixas com boa
rejeição de distúrbios e ruídos. No entanto, por questões de objetividade e
simplificação, não será abordada uma revisão sobre os aspectos probabilísticos do
Filtro de Kalman, onde estudam-se as características das perturbações e ruídos que
se deseja filtrar. Ao contrário, pretende-se ser mais direto ao se utilizar o
conhecimento sobre os distúrbios e ruídos mais comuns encontrados nos problemas
de controle de aeronaves. Para um estudo mais detalhado, recomenda-se a leitura
das referências Doyle e Stein (1979) e Stevens e Lewis (2003).
Quando não se tem conhecimento sobre as condições iniciais da aeronave
muito menos sobre as incertezas do modelo da planta e dos ruídos de medição, os
conhecimentos práticos da indústria aeronáutica podem fornecer informações que
preencham essas lacunas. Logo, de praxe, costuma-se projetar o Filtro de Kalman
para uma frequência de corte em torno de 10 rad/s para o máximo SV. Lembrando
também, que esta limitação de frequência adotada pelo filtro, algumas vezes entra
em conflito com o regulador LQR, fazendo com que o uso de ambos no LQG
necessite de alguma forma de recuperação, como o Loop-Transfer-Recovery.
4.5.4 LQG/LTR
Como foi visto nas seções anteriores deste capítulo, se um observador ou
Filtro de Kalman for utilizado para estimar os estados, as propriedades de garantia
de robustez do LQR com realimentação total de estados geralmente são perdidas
67
(Doyle e Stein; 1979). Desta forma, nesta seção deseja-se apresentar uma técnica
que recupera essa garantia de robustez e estabilidade. A técnica conhecida como
LTR, que recupera o loop-gain ou ganho de malha do regulador, pode ser feita a
partir da entrada ou da saída da planta. Neste trabalho, será abordada apenas a
técnica de recuperação do loop-gain na saída, que se baseia em projetar primeiro o
Filtro de Kalman e posteriormente o regulador LQR, ajustando os parâmetros deste
último e analisando as respostas no domínio do tempo e da frequência até que os
requisitos de robustez e estabilidade sejam atendidos.
A estrutura do controlador LQG/LTR pode ser vista na Figura 4.10. Esta
estrutura é a mesma a ser utilizada no projeto dos controladores laterais e
longitudinais do Aerosonde UAV neste trabalho.
Figura 4. 10: Diagrama de blocos do controlador LQG/LTR (Zarei et al.; 2006).
As matrizes Ac , Bc e Cc que aparecem no diagrama da Figura 4.10 correspondem às
matrizes Aa , Ba e Ca do controlador LQG/LTR descritas adiante na equação 4.32
para o caso lateral ou em (4.33) para o caso longitudinal. Outro bloco que deve ser
destacado é o do pré-compensador que aparece na Figura 4.10, pois é utilizado nos
projetos do Aerosonde UAV para balancear a dinâmica da planta, não sendo uma
parte integrante e necessária para o LQG/LTR.
O projeto do LQG/LTR para a recuperação da robustez do loop-gain G ( s) K ( s)
na saída, onde K ( s) é o controlador e G ( s) a planta, começa com o projeto do Filtro
de Kalman, obtendo assim o loop-gain de Kalman:
Lk ( s ) = C ΦL
(4.22)
onde Φ ( s ) = ( sI − A) −1 é a matriz de transição de estados. A equação 4.22 usufrui da
mesma propriedade de robustez garantida pela realimentação de estados do loop68
gain do LQR, K ΦB , quando a recuperação é feita na entrada. Isto se deve a
dualidade existente entre as técnicas de recuperação na entrada e na saída
(Stevens e Lewis; 2003).
Após o projeto do Filtro de Kalman, deve-se projetar o regulador LQR para
obter o loop-gain completo Lr ( s ) do regulador:
Lr ( s ) = G ( s ) K ( s ) = C ΦBK Φ r L
(4.23)
Φ r ( s ) = [ sI − ( A − BK − LC )]−1
(4.24)
onde Φ r ( s ) é a matriz de transição dos estados já com o Filtro de Kalman e o
regulador LQR, sendo necessário selecionar os parâmetros do LQR para que Lr ( s )
se aproxime de C ΦL . A função custo a ser minimizada pode ser da forma
J=
1 ∞ T
( x Qx + ρ 2uT Ru ) dt
2 ∫0
(4.25)
Após a seleção de Q e R , pode-se variar o valor de ρ , normalmente fazendo-o
tender a zero para efetuar o loop-shaping pela análise dos gráficos de magnitude
dos valores singulares e recalculando o ganho K até se obter resultado satisfatório.
O uso do LQG/LTR será empregado nas próximas seções para executar as
funções de controle de SAS, CAS e piloto automático do Aerosonde UAV. Mais
precisamente, na dinâmica lateral, serão projetados controladores que controlam o
ângulo de bank com aumento da estabilidade e promovendo um piloto automático do
tipo roll angle hold. Para a dinâmica longitudinal, o mesmo será feito para um pitch
angle hold com speed hold, todos em um único controlador com a estrutura sendo
fornecida automaticamente pelo processo, a seguir, do projeto dos controladores.
4.6 Sistema de Controle Lateral do Aerosonde UAV
O projeto do sistema de controle lateral do Aerosonde UAV começa pelo
desenvolvimento de um sistema CAS que controle o ângulo de bank da aeronave
para o problema de regulação e rastreamento, ou seja, tanto para manter as asas
niveladas como para mantê-la em ângulos específicos e permitir curvas
coordenadas com mínima derrapagem. Logo, o interesse está em fazer φ seguir um
comando do piloto automático mantendo β = 0 .
69
Este sistema será projetado utilizando as técnicas de LQG/LTR estudadas
anteriormente visando agrupar em uma única estrutura o roll-damper, yaw-damper,
roll-rate e o piloto automático roll-angle hold, que será o controlador de baixo nível
para o controle direcional do Aerosonde.
Utilizando o software Matlab para o projeto, e o modelo linear do Aerosonde,
apresentado no Capítulo 3, inicia-se esta etapa pela análise das variáveis que serão
utilizadas para estimar os estados da planta e comandar a mesma.
O vetor de entrada de controle é u = [ aileron rudder ] ;
T
O vetor de referência é r = ⎡⎣ rφ
O vetor de estados é x = [ v
T
rβ ⎤⎦ , onde rβ = 0 ;
p r φ ψ] ;
T
A matriz C da planta é:
⎡ 0
C=⎢
⎣0, 0435β
0 0 1φ
0 0
0
0⎤
0 ⎥⎦
logo, as variáveis que serão utilizadas para estimar os demais estados são φ e β ,
sendo feita a aproximação de β = 0, 0435v , obtida pelo software de linearização
descrito no ANEXO B, onde calcula-se o ângulo de sideslip a partir do vetor formado
pelas duas componentes u e v da velocidade nos eixos x e y, respectivamente.
Na Figura 4.11 é apresentado o diagrama de magnitude de Bode em relação
aos ganhos principais inferior ( σ ) e superior ( σ ) do sistema original da dinâmica
lateral do Aerosonde. Os ganhos principais (ou SVs) versus frequência mostram
claramente que o erro de regime será elevado em malha fechada, pois o loop-gain
não apresenta integradores nem um grande ganho DC, assim como uma grande
separação entre os SVs inferior e superior nas baixas freqüências. Logo, esta
característica deve ser balanceada pelo uso de um pré-compensador P :
P = H −1 (0)
(4.26)
que é a inversa do ganho DC H (0) do sistema original que pode ser obtido pelo
comando dcgain do Matlab:
⎡52,8559 48, 4368⎤
H (0) = ⎢
⎥
⎣ 1, 4301 2,3250 ⎦
e
⎡ 0, 0434 -0,9034 ⎤
P=⎢
⎥
⎣-0, 0267 0,9858 ⎦
70
O sistema aumentado pela inclusão do pré-compensador na entrada é descrito por
d ⎡ x⎤ ⎡ A B⎤ ⎡ x⎤ ⎡ 0 ⎤
=
+
u
dt ⎢⎣ε ⎥⎦ ⎢⎣ 0 0 ⎥⎦ ⎢⎣ε ⎥⎦ ⎢⎣ P ⎥⎦
Aa
;
⎡y⎤
⎡ x⎤
⎢ y ⎥ = [C 0] ⎢ε ⎥
Ca
⎣ ⎦
⎣ ε⎦
(4.27)
Ba
e o resultado em termos dos SVs versus frequência são apresentados na Figura
4.12. Percebe-se que o ganho DC dos SVs foi balanceado e agora o sistema é do
tipo 1, com integradores nas duas entradas, acarretando em uma inclinação de -20
dB/década, com ganho mais elevado nas baixas freqüências.
Margens de robustez no domínio da frequência
Considerando primeiro a margem de alta freqüência, assumindo que o
modelo do Aerosonde é preciso em 10% até a frequência de 2 rad/s. Após este valor
as incertezas aumentam com uma taxa de 20 dB/década devido aos erros de
modelagem, modos flexíveis da aeronave e ruídos de medição. Este comportamento
é modelado por:
m(ω ) =
s+2
20
(4.28)
Logo, para garantir a estabilidade da planta, o loop-gain referente a saída
deve ser:
σ (GK ( jω )) <
1
20
=
m(ω ) s + 2
(4.29)
mostrado na Figura 4.13. Sendo assim, o projeto do Filtro de Kalman será feito de
forma a respeitar a limitação da Figura 4.13, procurando atenuar comportamentos
acima de 10 rad/s.
Para as considerações de robustez em relação às baixas freqüências, o
sistema em malha fechada deve ser robusto, por exemplo, a distúrbios como rajadas
de ventos. O mínimo loop-gain para o SV inferior deve ser não muito menor que
20dB em 10−1 rad/s.
71
Figura 4. 11: Ganhos principais do sistema original da dinâmica lateral do Aerosonde.
Figura 4. 12: Ganhos principais balanceados em DC e com inclinação de -20 dB/década.
72
Figura 4. 13: Magnitude de
1
m (ω )
, para limitar
σ
nas altas freqüências.
4.6.1 LQG/LTR Lateral do Aerosonde
Para o projeto do LQG/LTR lateral com recuperação do loop-gain na saída,
primeiro será projetado o Filtro de Kalman, logo, C Φ ( s ) L é o que se quer recuperar
na fase posterior com o projeto do LQR.
Filtro de Kalman
Para encontrar o ganho L que minimiza o funcional
J=
1 ∞ 2 T
( rf x Qx + uT Ru ) dt
2 ∫0
(4.30)
os parâmetros de rf , Q e R devem ser selecionados para garantir rápida
convergência na estimação, boa resposta temporal e respeitar as margens de
robustez e estabilidade.
Para o sistema aumentado após a inclusão do pré-compensador P , de (4.26),
o vetor de estados aumentou seu número de variáveis para 7, logo a matriz de
pesos Q deve ser diagonal 7x7 e R , referente as duas variáveis de entrada, uma
73
matriz identidade 2x2. Após algumas tentativas, determinou-se Q do filtro como
sendo:
v
p
r
φ
ψ
P1
P2
Q = diag ([1 1 1 1 1 0,1 0,1])
e rf = 102 , o que forneceu os SVs da Figura 4.14, e resposta temporal da Figura 4.15
ao degrau unitário aplicado na referência rφ com o sistema em malha fechada.
Os valores singulares da Figura 4.14 estão respeitando as margens de
robustez e estabilidade, e a resposta temporal da Figura 4.15 é aceitável, com erro
de regime nulo.
O ganho L encontrado que minimiza o funcional em (4.30) é:
⎡ -6, 6814 17,1550 2,5413 11, 6021 1,5383 -2, 6822 -1, 6702 ⎤
L=⎢
⎥
⎣ 71, 2020 -7,5490 -6,8505 -0, 2906 -3, 4528 -1, 6704 2, 6829 ⎦
T
e as margens de ganho e fase do sistema são, respectivamente, 2 ×104 e 60º .
Outro requisito de análise para o projeto do filtro é a condição de
alcançabilidade de ( Aa , Q ) , que é atendida pela matriz de pesos Q . Logo, o
sistema é alcançável, o que pode ser verificado pelo erro de estimação tendendo a
zero na Figura 4.16 após a realização da simulação com o diagrama de blocos da
Figura 4.17.
Figura 4. 14: Ganhos principais de C Φ ( s ) L com Filtro de Kalman e rf = 10 .
2
74
Figura 4. 15: Resposta ao degrau aplicado em
rφ do sistema em M.F. com Filtro de Kalman.
Figura 4. 16: Erro de estimação dos estados tendendo a zero, exceto por ψ na linha preta tracejada.
75
Figura 4. 17: Diagrama de simulação do Simulink e análise do erro de estimação de estados.
Na Figura 4.16 a linha preta tracejada chama a atenção, pois é a única que
não converge totalmente a zero, deixando o erro de estimação em torno de -0,04.
Esta linha corresponde ao erro de estimação de ψ (psi). No entanto, esta variável
não causará problemas, pois ficará “relaxada” no projeto do LQR, justamente para
não influenciar no controle, já que está ligada diretamente ao controle de direção,
que não está sendo abordado neste controlador, para não entrar em conflito devido
ao acoplamento que ψ tem com todas as outras variáveis do vetor de estados.
No diagrama da Figura 4.17 destacam-se as matrizes A est , Best e Cest que são
as matrizes do Filtro de Kalman. Além destes destaques, a entrada aplicada aos
ailerons foi um pulso que durou de 0 a 1 segundo de simulação com amplitude igual
a 0,1 apenas para tirar a planta do seu estado inicial de equilíbrio.
Os pólos, fatores de amortecimento, freqüências ou constantes de tempo com
o filtro em malha fechada são mostrados na Tabela 4.2 (veja especificações de
qualidade na seção 4.1).
Tabela 4. 2: Pólos de malha fechada do sistema com o Filtro de Kalman.
Pólos
Amortecimento
Freq. (rad/s) ou Cte. Tempo (s)
-1,53
-2, 07
-1,94 ± j 5, 60
-9,53
-18,3
1, 00
1, 00
0,328
1, 00
1, 00
1,53 s
2, 07 s
5,93 rad / s
9,53 s
18,3 s
76
LQR
A seleção dos parâmetros do LQR, que leva ao ganho ótimo K , que minimiza
o funcional
J=
(
)
1 ∞ T
x Qx + ρ 2u T Ru dt
2 ∫0
(4.31)
deve ser realizada. O parâmetro ρ 2 = 10−2 foi obtido após algumas tentativas e
análises da resposta do sistema com o LQG/LTR completo. Escolheu-se a matriz de
pesos Q do LQR dada por:
φ
ψ
P1
P2
⎛ v 2 p r
⎞
6
−6
Q = diag ⎜ ⎡⎣10 1 1 10 10
1 1⎤⎦ ⎟
⎜
⎟
⎝
⎠
Percebe-se que esta matriz Q deixa a variável ψ bastante “relaxada” tal como
mencionado anteriormente sobre a Figura 4.16. No entanto, a variável v , que tem
relação direta com a variação de β , possui maior peso assim como a variável φ . A
matriz R , é uma matriz identidade 2x2, agora multiplicada por ρ 2 = 10−2 . Os testes e
análises a partir desta etapa já são feitos com o sistema completo do LQG/LTR, e o
valor calculado é dado por:
⎡0, 0849 -0, 0849 -0, 4828 -3, 4213 0, 0000 1, 0272 0,9352 ⎤
K = 103 ⎢
⎥
⎣ 0, 0271 0,1539 -0, 0604 9,3653 0, 0000 -0, 0705 0,1137 ⎦
LQG/LTR
O projeto do LQR mostrado anteriormente já atende as necessidades de
recuperação do loop-gain, ou seja, a variação do parâmetro ρ , assim como das
matrizes dos pesos, já foram feitas para garantir a recuperação das especificações
de C ΦL . Sendo então apresentado o sistema aumentado e completo com o
LQG/LTR, pré-compensador e o modelo linear do Aerosonde:
⎡A
Alqg / ltr = ⎢ a
⎣ 07 x 7
− Ba K
⎤
⎥
Aa − LCa − Ba K ⎦
;
⎡0 ⎤
Blqg / ltr = ⎢ 7 x 2 ⎥
⎣ −L ⎦
;
Clqg / ltr = [Ca
02 x 7 ] (4.32)
Os valores singulares do sistema completo, respeitando as margens de robustez e
estabilidade, são mostrados na Figura 4.18, e a resposta a aplicação de um degrau
unitário na referência rφ e o sinal de controle, na Figura 4.19.
77
As margens de ganho e fase do sistema com LQG/LTR são respectivamente
6, 2001×104 e 60, 0011º , fornecendo um aumento na margem de ganho em relação a
obtida no projeto com o Filtro de Kalman.
Figura 4. 18: Ganhos principais do sistema com LQG/LTR.
Figura 4. 19: Resposta ao degrau do sistema com LQG/LTR e sinal de controle.
Na Figura 4.19 vê-se que o sinal de controle conta com baixas amplitudes
angulares, levando em consideração que amplitudes comuns de deflexão dos
ailerons e do leme podem chegar a até 45º em algumas aeronaves.
78
O controlador lateral de baixo nível está pronto, então, para receber os
comandos do controlador de direção que vai orientar o Aerosonde em relação aos
ângulos da bússola através de comandos de ângulos de rolagem ao piloto
automático Roll Angle Hold (controle do ângulo de bank) projetado aqui.
Como uma observação final, percebe-se, na Figura 4.19, que a resposta
temporal do ângulo de rolagem possui sobre-sinal de aproximadamente 10% que
não foi possível reduzir sem comprometer as margens de robustez e estabilidade.
No entanto, em aviões, os efeitos que mais provocam desconforto a tripulantes ou
danos aos equipamentos dentro da fuselagem são os de derrapagem e demais
esforços resultantes das forças que provocam o efeito de sideslip (Stevens e Lewis;
2003). Portanto, sobre-sinal nesta faixa para o Roll Angle Hold é aceitável, tal como
pode ser comparado pela Tabela 4.3 dos pólos de malha fechada do sistema, em
relação às especificações de qualidade de vôo da seção 4.1 e os pólos de malha
aberta do sistema original na Tabela 3.2 do Capítulo 3.
Tabela 4. 3: Em cinza, os pólos de M.F. do Aerosonde aumentados pelo pré-compensador e modificados pelo
controlador. Em branco, os pólos adicionados pelo LQG/LTR.
Pólos
Amortecimento
Freq. (rad/s) ou Cte. Tempo (s)
−1,86 × 10−6
1
1,86 ×10−6 s
−1,53
1
1,53 s
−2, 07
1
2, 07 s
−1,94 ± j 5, 60
0,328
5,93 rad / s
−9,53
1
9,53 s
−4,84 ± j8,31
0,504
9, 62 rad / s
−9, 66
1
9, 66 s
−18,3
1
18,3 s
−49,3 ± j83, 7
0,507
97, 2 rad / s
−98, 0
1
98, 0 s
Todo o código desenvolvido no Matlab para o projeto do Controlador Lateral
do Aerosonde, que foi apresentado nesta seção, está presente no APÊNDICE B.
79
4.7 Sistema de Controle Longitudinal do Aerosonde UAV
O projeto do sistema de controle longitudinal do Aerosonde UAV incorpora os
pilotos automáticos de Pitch Angle Hold e Speed Hold para o controle das variáveis
θ e Va (velocidade do fluxo de ar), respectivamente. O método de projeto se
assemelha ao utilizado para o controlador LQG/LTR lateral de baixo nível da seção
4.6, logo, este sistema incorpora aumento de estabilidade com SAS, CAS e funções
de piloto automático.
No controle longitudinal, uma das maiores preocupações é o desacoplamento
das variáveis que se quer controlar, pois como foi visto no Capítulo 3, velocidade e
atitude da aeronave estão fortemente ligadas. Além desta preocupação, há também
a questão de minimizar as variações angulares e forças laterais que atuam na
dinâmica longitudinal e “relaxar” algumas variáveis de estados, tal como a altitude e
a rotação do motor (Omega), para que estas não interfiram no controle de θ e Va .
O vetor de entrada de controle é: u = [ elevator throttle] ;
T
O vetor de referencia é: r = [ rθ
rVa ] ;
T
O vetor de estados é: x = [u w q θ
h omega ] , onde os componentes
T
são, respectivamente: velocidade ( m / s ) em relação ao eixo x, velocidade ( m / s ) em
relação ao eixo z (sistema NED), pitch-rate, ângulo de pitch, altitude e Omega (rpm);
A matriz C da planta é:
⎡0,9972Va
C=⎢
0
⎣
0, 0753α
0
0 0
0 1θ
0 0⎤
0 0 ⎥⎦
logo, as variáveis que serão utilizadas para estimar os demais estados são
Va , α e θ .
Na Figura 4.20 é apresentado o diagrama de magnitude de Bode em relação
aos ganhos principais inferior ( σ ) e superior ( σ ) do sistema original da dinâmica
longitudinal do Aerosonde e, assim como ocorreu com a dinâmica lateral, será
necessário incluir um pré-compensador P como o da equação 4.26, para a inclusão
dos integradores e balancear as entradas do sistema. Neste caso, P é dado por:
⎡-0, 0006 -3, 6743 ⎤
P=⎢
⎥
⎣ 0,1929 12,1762 ⎦
80
O sistema aumentado com P é descrito da mesma forma que em (4.27).
Figura 4. 20: Ganhos principais inferior e superior do sistema longitudinal original.
Figura 4. 21: Ganhos principais do sistema aumentado com pré-compensador.
81
Na Figura 4.21, são mostrados os valores singulares do sistema aumentado
com o pré-compensador. Percebe-se que agora o sistema decai com -20 dB/década,
mas diferentemente do que ocorreu com a dinâmica lateral, os canais não estão
perfeitamente balanceados. No entanto, houve um aumento significativo do ganho
em DC.
4.7.1 LQG/LTR Longitudinal do Aerosonde
Filtro de Kalman
Para encontrar o ganho L que minimiza o funcional em (4.30) foram
utilizados, neste caso, os seguintes parâmetros:
rf 2 = 10
u
w
q
θ
h omega P1
P2
Q = diag ([1 1 1 1 1 1 30 90])
R = I2x2
resultando em:
⎡ 0, 0111 -0, 0238 -0, 0041 -0, 0024 0, 0087 2, 7802 0, 0036 0, 0293⎤
L = 10 ⎢
⎥
⎣-0, 0096 0, 0956 0, 0239 0, 0072 0, 0028 0, 6060 -0, 0169 0, 0062 ⎦
T
3
Sendo o sistema alcançável e observável pelos testes de ( Aa , Q ) e ( AaT , BaT ) .
Os SVs do sistema em malha aberta com o Filtro de Kalman são mostrados
na Figura 4.22, projetados para respeitar os limites de altas freqüências tal como foi
feito com o controlador lateral, tendo frequência de corte em torno de 10 rad / s .
Para avaliar o funcionamento do filtro, na Figura 4.23 é apresentada a
resposta temporal quando um degrau é aplicado à referência rθ , com rVa = 0 . Na
Figura 4.24 é feito o contrário, aplicando-se o degrau em rVa e zero na referência 1.
O erro de estimação é avaliado na Figura 4.25, quando um diagrama de simulação
tal como o apresentado na seção anterior, na Figura 4.17, é utilizado, aplicando-se
um pulso de amplitude igual a 0,1 na entrada (1) do elevator durante o primeiro
segundo da simulação.
As margens de ganho e fase do sistema com o filtro são respectivamente
1,9804 × 105 e 60º .
82
Figura 4. 22: Ganhos principais do sistema aumentado em M.A. com o Filtro de Kalman.
Figura 4. 23: Resposta temporal ao degrau aplicado em rθ do sistema em M.F. com o Filtro de Kalman.
83
Figura 4. 24: Resposta temporal ao degrau aplicado em rV a do sistema em M.F. com o Filtro de Kalman.
Figura 4. 25: Erros de estimação dos estados do sistema com o Filtro de Kalman.
84
Na Figura 4.25, pode-se perceber que uma das variáveis de estados fica com
erro estático não nulo. Essa variável é referente à altitude, sobrando um erro de
estimação de aproximadamente 0, 244 metros. Como mencionado na seção 4.4.2,
deste capítulo, as tolerâncias nos controladores de altitude são de ±60 m , logo, este
erro de estimação não é tão significativo. Além disto, a variável h , assim como
Omega , ficarão “relaxadas” no projeto deste controlador, justamente para não
interferirem no objetivo de controle deste LQG/LTR longitudinal.
Os pólos, fatores de amortecimento, freqüências e constantes de tempo com
o filtro em malha fechada são mostrados na Tabela 4.4 (veja especificações de
qualidade na seção 4.1).
Tabela 4. 4: Pólos de M.F. do sistema aumentado com o Filtro de Kalman.
Pólos
Amortecimento
Freq. (rad/s) e Cte. Tempo (s)
-2,89 ×10−4
-2, 67 ± j 0,991
1
0,937
2,89 × 10−4 s
2,85 rad / s
-4,89
-4, 42 ± j 3,97
-4, 28 ± j10,3
1
0, 744
4,89 s
5,94 rad / s
11,1 rad / s
0,385
LQR
Para a obtenção do ganho K ótimo que minimiza o funcional em (4.31), os
seguintes parâmetros foram utilizados no projeto do LQR:
ρ 2 = 10−2
q
θ
h
omega
P1
P2
⎛ u 3 w
⎞
5
−3
−3
−3
−3
⎡
Q = diag ⎜ ⎣5 ×10 10 1 10 10
10
10
10 ⎤⎦ ⎟
⎜
⎟
⎝
⎠
R = I2x2
Percebe-se que a matriz Q deixa as variáveis h e Omega “relaxadas”, enquanto que
u , que tem relação direta com Va , possui um peso bem maior. A variável θ também
possui grande peso, assim como w e q também recebem certa consideração para
minimizar as respectivas variações de velocidade no eixo z e variação angular de θ .
O ganho K obtido para estes parâmetros de projeto foi:
85
⎡ 0, 6580 0, 0909 -0, 0545 -1,3108 0, 0000 0, 0013 0, 2307 0, 0696 ⎤
K = 103 ⎢
⎥
⎣-0,1031 0, 0004 0, 0873 3, 0785 0, 0003 0, 0001 -0, 0338 0, 0009 ⎦
LQG/LTR
A descrição do sistema aumentado e completo com o LQG/LTR, précompensador e o modelo linear longitudinal do Aerosonde é:
⎡A
Alqg / ltr = ⎢ a
⎣08 x8
− Ba K
⎤
Aa − LCa − Ba K ⎥⎦
;
⎡0 ⎤
Blqg / ltr = ⎢ 8 x 2 ⎥
⎣ −L ⎦
;
Clqg / ltr = [Ca
02 x8 ]
(4.33)
Os valores singulares do sistema completo, respeitando as margens de robustez e
estabilidade são mostrados na Figura 4.26, e a resposta à aplicação de um degrau
unitário na referência rθ e o respectivo sinal de controle, na Figura 4.27. A resposta
à aplicação de um degrau unitário na referência rVa e o respectivo sinal de controle
são mostrados na Figura 4.28.
As margens de ganho e fase do sistema com LQG/LTR são respectivamente
1,1645 ×105 e 60º .
Figura 4. 26: Ganhos principais do sistema com LQG/LTR.
86
Na Figura 4.26 vê-se que os ganhos principais inferior e superior respeitam as
margens de robustez e estabilidade estabelecidas para os projetos dos
controladores LQG/LTR, tendo o ganho principal inferior em 0,1 rad / s em torno de
20dB e o superior atenuando a partir de 10 rad / s .
Figura 4. 27: Resposta ao degrau aplicado em rθ do sistema em M.F. com LQG/LTR e o sinal de controle.
Figura 4. 28: Resposta ao degrau aplicado em rVa do sistema em M.F. com LQG/LTR e o sinal de controle.
A resposta de θ na Figura 4.27 possui um elevado sobre-sinal de cerca de
60%, considerado ruim para grandes deflexões do ângulo de pitch. No entanto, o
acoplamento existente entre Va e θ no Aerosonde UAV não deixou muita escolha
para garantir a robustez e estabilidade de acordo com as margens estabelecidas no
projeto. De qualquer forma, será visto na simulação com a planta não-linear, que
com a inclusão de um limitador de pitch (no Bloco 3 da Figura 5.7), este controlador
87
pode ser limitado a rastrear apenas valores pequenos de ângulo de pitch. Além
disto, será visto no Capítulo 5 que os resultados com o modelo não linear possuem
menor sobre-sinal. Isto se deve ao fato de que a dinâmica dos atuadores não é
considerada no modelo linear, assim como efeitos de saturação. Este fato pode ser
explicado pelos gráficos dos sinais de controle das figuras 4.27 e 4.28, onde
percebe-se variações muito rápidas para as deflexões angulares do profundor. O
mesmo ocorre para o acelerador, que no modelo linear, nunca satura, o que permite
imaginar uma alavanca de acelerador infinita que aceita até mesmo valores
negativos (desaceleração ou freio). No modelo não linear o acelerador vai de 0 a 1.
Por hora, apesar do sobre-sinal estar elevado, sabe-se que o sistema é
robusto e estável, sendo capaz de garantir um bom desacoplamento entre a
dinâmica longitudinal e lateral do Aerosonde, garantindo erro de regime nulo durante
curvas coordenadas quando trabalhando em conjunto com o controlador lateral.
Os pólos do sistema longitudinal completo são mostrados na Tabela 4.5:
Tabela 4. 5: Em cinza, os pólos M.F. do Aerosonde aumentados pelo pré-compensador e modificados pelo
controlador. Em branco, os pólos adicionados pelo LQG/LTR.
Pólos
Amortecimento
Freq. (rad/s) ou Cte. Tempo (s)
-2,89 ×10−4
1, 00
2,89 × 10−4 s
-7,90 × 10−4
1, 00
7,90 ×10−4
-2, 67 ± j 0,991
0,937
2,85 rad / s
−3,53
1
3,53 s
−4,89
1
4,89 s
-4, 42 ± j 3.97
0, 744
5,94 rad / s
−7, 70
1
7, 70 s
-4, 08 ± j 6, 61
0,525
7, 77 rad / s
-4, 28 ± j10,3
0,385
11,1 rad / s
-99,8
1
99,8 s
-50,8 ± j87, 4
0,502
101 rad / s
O código em Matlab utilizado neste projeto está descrito no APÊNDICE C.
O projeto dos sub-controladores do Aerosonde está finalizado, passando-se
então ao projeto dos sistemas de controle de altitude, direção e navegação.
88
Capítulo 5
Controladores fuzzy de Altitude, de Direção e o Sistema de Navegação
No Capítulo 4 foram desenvolvidos os controladores de baixo nível do sistema
de controle de vôo do Aerosonde UAV, sendo estes, sub-controladores, que são
utilizados pelos controladores fuzzy de altitude e de direção apresentados neste
capítulo. Além desses, será mostrado o sistema de navegação da aeronave que
fornece as referências aos controladores de altitude e direção para coordenar o UAV
através dos waypoints pré-estabelecidos.
A escolha em trabalhar com controle fuzzy, nesta etapa, baseia-se em
questões didáticas de apreciação do uso concomitante deste tipo de técnica com
outras, e pela simplicidade na elaboração dos controladores, em vista de que os
acoplamentos e demais complexidades da planta já foram solucionadas pelo método
LQG/LTR dos controladores do Capítulo 4. Além disto, este trabalho complementa o
que foi desenvolvido sobre controle fuzzy aplicado ao Aerosonde em Silveira et al.
(2007), onde não se fazia o controle simultâneo da altitude e velocidade, nem a
navegação através dos waypoints.
Ao final do capítulo, serão feitas simulações com o sistema de controle de vôo
completo juntamente com o sistema de navegação. Serão verificados resultados
com e sem perturbações, testes fora da região de estabilidade onde o modelo
utilizado nos projetos dos controladores foi obtido por linearização, e também com
variações na massa de combustível. Além disso, serão fornecidas informações sobre
cada bloco do Simulink desenvolvido para simular este sistema, com o modelo nãolinear com 6-DOF.
5.1 Controlador fuzzy de Altitude
As técnicas de controle fuzzy estão inseridas na área de controle inteligente.
Geralmente, os problemas de controle são tratados através destas técnicas quando
um modelo qualitativo está disponível (ou informações de especialistas) e a
estratégia de controle pode ser formulada e executada com base em um conjunto de
regras que procuram emular o raciocínio humano.
89
Uma das principais características dos controladores fuzzy do tipo Mamdani é
a presença do fuzzificador e defuzzificador que permitem trabalhar com regras
lingüísticas, tal como, “SE” a altitude é grande, “ENTÃO” diminua o ângulo de pitch.
Essa expressão SE-ENTÃO mostrada caracteriza uma regra que qualquer piloto de
avião conhece e faz intuitivamente. Sendo assim, o que se pretende fazer para
controlar a altitude do Aerosonde UAV é utilizar uma base de conhecimento para
promover o controle.
Como os sinais fornecidos pelos sensores do Aerosonde são valores reais, ou
crisp (Passino e Yurkovich; 1998), deve-se utilizar um fuzzificador na entrada do
controlador para transformar um valor real em um conjunto fuzzy. O processo
inverso é feito na saída com o defuzzificador, transformando um conjunto fuzzy de
volta para um valor real, como por exemplo, um ângulo de pitch. A estrutura em
diagrama de blocos de um controlador fuzzy de Mamdani é mostrada na Figura 5.1
(Wang; 1997):
Figura 5. 1: Diagrama de blocos do controlador fuzzy de Mamdani.
Para melhor entender o processo mostrado na Figura 5.1, a elaboração do
piloto automático Altitude Hold, abordado na seção 4.4.2, será usado como exemplo.
O controlador fuzzy de altitude trabalho com mapeamento estático, recebendo
na entrada o erro entre a altitude desejada e a atual do Aerosonde UAV
eh = href − h
(5.1)
para determinar uma referência de ângulo de pitch ( θ ref ) na saída para fazer a
aeronave subir, nivelar ou descer. O processo de fuzzificação então converte o valor
real de eh para os conjuntos fuzzy N, Z, P, referentes a erro negativo, zero e positivo,
respectivamente. Estes três conjuntos são representados por funções de pertinência
definidas no universo de discurso de eh , conforme, por exemplo, a Figura 5.2. A
fuzzificação é feita, então, mapeando-se o valor real da variável eh em termos das
funções trapezoidais de N e P ou da função triangular de Z (Wang; 1997), tal como
pode ser visto no código fonte do controlador fuzzy de altitude no APÊNDICE D.
90
Para o processo de defuzzificação, é aplicado o método por média dos
centros e são utilizados os conjuntos fuzzy referentes aos ângulos de pitch que
devem conduzir a aeronave a ganhar ou perder altitude. Estes conjuntos são
designados por AUMENTAR, NEUTRALIZAR e DIMINUIR o ângulo de pitch, com
suas funções, neste trabalho, mostradas na Figura 5.3.
Figura 5. 2: Funções de pertinência de entrada do controlador.
Figura 5. 3: Funções de pertinência de saída do controlador.
As funções de pertinência da saída da Figura 5.3 estão limitadas em uma
faixa de 0 a 1 no universo de discurso de θ ref , no entanto, o valor crisp resultante da
defuzzificação é posteriormente escalonado por um ganho que limita o θ ref máximo
para que se possa definir a máxima deflexão de pitch de acordo com a fase de vôo,
como por exemplo, para não permitir grandes ângulos em baixas velocidades. A
distribuição das funções de pertinência tanto da entrada como da saída são
geralmente determinadas com base no conhecimento adquirido sobre a planta.
91
A máquina de inferência aqui utilizada usa uma base de regras simples
composta por:
SE eh ∈ N ENTÃO diminuir pitch;
SE eh ∈ Z ENTÃO neutralizar pitch;
SE eh ∈ P ENTÃO aumentar pitch.
Como se usou apenas uma entrada, a inferência é direta tal como mostrado na base
de regras. E para a defuzzificação, utiliza-se a média dos centros:
3
θ ref =
∑b μ
i =1
3
i
∑μ
i =1
i
(5.2)
i
onde b é a base da função de saída. Neste exemplo, pela análise de Figura 5.3 vêse que as bases são 0, 25 , 0,5 e 0, 75 para diminuir, neutralizar e aumentar,
respectivamente, e μ é a pertinência das funções de entrada. O código fonte deste
controlador está demonstrado no APÊNDICE D.
5.2 Controlador fuzzy de direção/heading
O controlador de direção ou Heading Hold é a parte do piloto automático
responsável por rastrear valores de ângulos entre 0º e 359º fornecidos pelo sistema
de navegação do Aerosonde. Na sua entrada recebe o erro entre a direção desejada
e a atual do UAV:
eψ = ψ ref −ψ UAV
(5.3)
ψ UAV = ψ + β
(5.4)
onde
A equação 5.4 representa a direção real da aeronave considerando os efeitos de
derrapagem ( β ). Como o controlador lateral desenvolvido no Capítulo 4 elimina
quase que totalmente os efeitos de β , pode-se considerar este igual a zero em
(5.4).
Na saída do controlador fuzzy de direção é fornecido o ângulo de bank ( φref )
de referência para o controlador lateral de baixo nível.
92
A estrutura do controlador fuzzy de direção é a mesma utilizada para o
controlador de altitude da seção anterior, diferenciando-se é claro, pelas variáveis de
entrada e saída assim como pela distribuição das funções de pertinência nos seus
universos de discurso. Para a entrada, os conjuntos continuam com os mesmos
nomes, N, Z e P. Para a saída, os novos nomes são: ESQUERDA, NIVELAR e
DIREITA, referentes a efetuar rolagem à esquerda ou à direita para mudar a direção
do UAV ou nivelar, para manter o veículo na direção atual. As funções de pertinência
para a entrada e saída são mostradas nas figuras 5.4 e 5.5.
Figura 5. 4: Funções de pertinência de entrada do controlador.
Figura 5. 5: Funções de pertinência da saída do controlador.
A saída do controlador fuzzy é o ângulo de bank de referência que é passado
ao controlador de direção de baixo nível, sendo que o universo de discurso de saída
está limitado entre -1 e 1, para posteriormente ser escalonado por um ganho externo
que limita o máximo ângulo de bank.
A base de regras do controlador fuzzy de direção é:
SE eψ ∈ N ENTÃO rolar à esquerda;
93
SE eψ ∈ Z ENTÃO nivelar asas;
SE eψ ∈ P ENTÃO rolar à direita.
O método de defuzzificação utilizado é o mesmo do controlador de altitude,
por média dos centros, mostrado na equação 5.2.
Deve-se fazer uma observação sobre o controlador de direção, pois ele
possui um tratamento especial para o sinal do erro de direção ( eψ ), para que a
aeronave sempre procure o caminho mais curto no momento de decidir para qual
lado efetuar uma curva. Por exemplo, imagine que o UAV esteja voando na direção
de ψ UAV = 10º quando o sistema de navegação solicita uma nova referência de
ψ ref = 350º . A aeronave poderia alcançar o ângulo solicitado tanto com uma curva à
esquerda como à direita, no entanto, pela esquerda o caminho é mais curto. Para
tratar este problema foi utilizado o esquema da equação 5.5:
(
)
eψ corrigido = eψ − 180º − 180º
se
eψ > 180º
(5.5)
A implementação de (5.5) assim como do controlador fuzzy de direção pode ser
conferida no APÊNDICE E.
Com a conclusão do controlador de direção, completa-se todo o sistema de
controle de vôo do Aerosonde UAV, restando apenas o sistema de navegação que
monitora a aeronave e comanda os controladores desenvolvidos neste trabalho.
5.3 Sistema de navegação
O sistema de navegação implementado neste trabalho deve exercer as
funções tal como abordado na seção 2.3.3. Basicamente, monitora-se a posição da
aeronave em termos de sua latitude, longitude e altitude determinando a partir de
uma base de dados de waypoints como conduzir o avião ao ponto desejado.
Neste trabalho, os waypoints possuem informações adicionais para que se
possa ter condições de vôo específicas entre cada fase da viagem. Essas
informações são velocidade, pitch máximo e bank máximo, permitindo que
decolagens, subidas até altitude de cruzeiro, descidas, aproximações e pousos
possam ter configurações que garantam a segurança da aeronave. Um exemplo da
utilização destes parâmetros adicionais dos waypoints seria ter o Aerosonde, logo
94
após a decolagem, sujeito a fortes ventos laterais e executando manobras com
grandes deflexões angulares. Nesta situação, se o avião ainda estivesse em baixa
velocidade e altitude, possivelmente poderia entrar em stall e cair.
Apesar dos parâmetros adicionais presentes na base de dados de waypoints
do sistema de navegação, este trabalha com restrições apenas de latitude e
longitude, ou seja, se a latitude e longitude de um determinado ponto for atingida, o
sistema de navegação passa a conduzir a aeronave para o próximo ponto desejado,
independentemente dos demais parâmetros estarem ou não próximos aos valores
desejados. Isto foi adotado por questões de simplificação do trabalho.
O método para determinação da trajetória que conduz o Aerosonde para uma
posição desejada em termos de latitude e longitude é o de triangulação, semelhante
ao adotado em Hoover (2004), onde
⎛ Lond − Lonatual ⎞ 180º
⎟
⎝ Latd − Latatual ⎠ π
ψ ref = tan −1 ⎜
(5.6)
é obtido a partir da Figura 5.6. O subscrito “d” em (5.6) é referente à latitude e
longitude desejadas. O ângulo ψ ref (em graus) obtido é então passado como
referência a ser rastreada pelo controlador fuzzy de direção.
Figura 5. 6: Representação gráfica do método de triangulação para obtenção da trajetória ψ ref .
95
O sistema de navegação, então, utiliza duas entradas (latitude e longitude)
que são constantemente avaliadas para poder determinar a trajetória para um
determinado waypoint. Estas entradas são fornecidas como ângulos em radianos e
comparadas com a latitude e longitude desejadas. Quando
| Latd − Latatual | < 10−5 rad e
| Lond − Lonatual | < 10−5 rad
(5.7)
o sistema de navegação considera que o waypoint fora atingido e deve-se passar ao
waypoint seguinte.
Na saída do sistema de navegação são fornecidas as referências para os
sistemas de controle de vôo conforme mostrado no vetor em (5.8):
saídanav = ⎡⎣ψ ref
href − 3 Varef
φMAX
θ MAX
wp ⎤⎦
T
(5.8)
onde wp é o número do waypoint para onde o sistema de navegação estará
conduzindo o Aerosonde. Percebe-se que a referência para a altitude href
é
subtraída de 3 metros. Isto se deve a necessidade de corrigir o erro em regime
permanente que o controlador fuzzy de altitude não elimina por não possuir ação
integral, sendo então, mais simples, incluir este offset de -3.
Para o planejamento das missões que se quer executar com o Aerosonde e
definição dos waypoints, pode-se utilizar coordenadas para longitude e latitude
definidas em graus, minutos e segundos convertendo-as para radianos usando a
função dms2rad do Matlab, tal como é descrito no código fonte do planejador de
missões do APÊNDICE F. As coordenadas podem ser obtidas a partir de um GPS,
mapa, cartas de navegação, etc. Neste trabalho, foram utilizados os softwares
Google Earth (http://earth.google.com/) e o Microsoft Flight Simulator 2004 para
obter coordenadas sobre pontos específicos do globo terrestre.
O código fonte do sistema de navegação é apresentado no APÊNDICE G
deste trabalho.
5.4 Simulações com o piloto automático completo
Neste estágio, tudo o que era necessário para controlar e navegar o
Aerosonde de forma autônoma já foi projetado e está pronto para ser utilizado com o
modelo não-linear de 6-DOF da aeronave no Simulink.
96
As simulações serão feitas em etapas que seguem desde o planejamento da
missão, testes individuais dos controladores, até o momento em que todos são
interligados e a missão é executada do início ao fim com a verificação da rota
percorrida. No total, duas simulações serão executadas, sendo a primeira, com o
Aerosonde nas proximidades da região de estabilidade em que os controladores de
baixo nível foram projetados. Na segunda, o Aerosonde sofrerá diversas variações
de parâmetros e condições de vôo, sendo também submetido a perturbações como
rajadas de ventos. Ao final, os resultados dessas duas simulações serão
comparados para avaliar a eficácia dos controladores desenvolvidos.
5.4.1 Simulação 1: na região de estabilidade e sem perturbações
Primeiramente, deve-se analisar o diagrama de simulações construído no
Simulink e seus sub-blocos para compreender e associar suas funções com os
sistemas que foram projetados. Na Figura 5.7 é apresentado o diagrama completo,
enumerado para facilitar a explicação, conforme:
Bloco 1: Controlador lateral do Aerosonde desenvolvido na seção 4.6 para
eliminar os efeitos de derrapagem ( β → 0 ) e rastrear um ângulo de bank ( φref ). O
diagrama contido dentro do bloco 1, com a estrutura do LQG/LTR, é mostrado na
Figura 5.8. Os sinais de controle gerados pelo controlador devem ser convertidos de
graus para radianos antes de serem inseridos no modelo não-linear do Aerosonde;
Bloco 2: Controlador longitudinal do Aerosonde desenvolvido na seção 4.7
para controlar a velocidade do fluxo de ar ( Va ) e o ângulo de pitch ( θ ). O diagrama
interno do bloco 2 é apresentado na Figura 5.9. Atenção especial deve ser dada
para a presença da saturação de [0,1] e ganho de 1/ 349, 2659 que limitam o sinal de
comando enviado ao throttle. O ganho elimina a dinâmica de Omega (rotação do
motor) considerado na matriz B do modelo linear da seção 3.3.1.1. Além disto, o
sinal de controle enviado à superfície aerodinâmica sofre conversão para radianos;
97
Bloco 3: Controlador fuzzy de altitude da seção 5.1, com seus sub-blocos
mostrados na Figura 5.10. Utiliza o bloco Matlab Function do Simulink para carregar
a função FLC_altitude.m que tem o código fonte descrito no APÊNDICE D. Na saída
da função é utilizado um ganho de 10 / 7, 491 para escalonar o sinal de controle para
que seja possível selecionar um θ MAX ;
Bloco 4: Controlador fuzzy de direção desenvolvido na seção 5.2. Seus sub-
blocos são semelhantes aos do bloco 3, carregando neste a função FLC_heading.m
mostrada no APÊNDICE E. Os sub-blocos são apresentados na Figura 5.11 e
destaca-se o bloco Pi Bound que converte os sinais de entrada que chegam entre 0º
e 359º para [−180º ,180º ] . Na saída da função FLC_heading é utilizado um ganho de
1/ 0,5 para escalonar o sinal de referência do ângulo de bank que é passado ao
bloco 2 limitado pelo φMAX ;
Bloco 5: Sistema de navegação com seus sub-blocos apresentados na
Figura 5.12. Trata os sinais dos sensores do Aerosonde utilizando demultiplexadores
para obter apenas as informações de posição do GPS da aeronave. Para a obtenção
da trajetória e gerenciamento dos waypoints o bloco Matlab Function é carregado
com o psi_waypoint.m descrito no APÊNDICE G. Duas saídas da função
psi_waypoint.m não aparecem no diagrama, que são θ MAX e φMAX . Elas são variáveis
globais que são atualizadas internamente no Matlab durante a simulação;
Bloco 6: Bloco que acompanha a versão 1.2 do Aerosim Blockset para
interligar o Simulink ao Microsoft Fligth Simulator. Este bloco passa os parâmetros
de posição, ângulos de Euler e Va via o protocolo do software FSUIPC desenvolvido
por Peter Dowson (http://www.schiratti.com/dowson.html).
Bloco 7: O gerenciador de saídas que apenas armazena as diversas saídas
do Aerosonde em variáveis globais do Matlab para que seja reduzido o número de
linhas de sinais no diagrama de simulações. Para isto, foram usados os blocos
From-To do Simulink. As variáveis deste bloco são mostradas na Figura 5.13.
98
Figura 5. 7: Diagrama completo de simulação no Simulink.
99
Figura 5. 8: Diagrama interno do bloco 1.
Figura 5. 9: Diagrama interno do bloco 2.
Figura 5. 10: Diagrama interno do bloco 3.
Figura 5. 11: Diagrama interno do bloco 4.
100
Figura 5. 12: Diagrama interno do bloco 5.
Figura 5. 13: Diagrama interno do bloco 7.
5.4.1.1 Planejamento da missão
Utilizando o software Google Earth foram obtidas coordenadas para efetuar
uma missão sobre a cidade de Belém-PA, Brasil, que será chamada de Sobrevôo
em Belém. O Aerosonde deverá passar por cinco waypoints iniciando sua viagem
sobre o Aeroporto Júlio César em Belém a 1000 metros de altitude. A rota que a
aeronave deverá percorrer é apresentada na Figura 5.14.
101
Figura 5. 14: Rota que deverá ser percorrida passando por cinco waypoints.
O Aerosonde vai iniciar sua viagem na direção ψ 330º , alinhado com a
cabeceira da pista 33 do Aeroporto Júlio César na altitude, velocidade e massa de
projeto, tal como o modelo linear foi obtido. Seguindo na direção de 330º ,
interceptará o WP 1 que fica sobre a outra extremidade da pista do aeroporto.
Efetuando uma curva coordenada à esquerda, seguirá até a Baía do Guajará
interceptando o WP 2. No caminho até o WP 3, o Aerosonde passará nas
proximidades das Docas do Estado do Pará, Mercado do Ver-o-Peso e outros
pontos turísticos da cidade de Belém do Pará até iniciar seu retorno ao aeroporto via
os WP 4 e 5. A seguir, são apresentados os waypoints carregados para o sistema de
navegação e plotados na Figura 5.15 pelo software waypoint_gen.m descrito no
APÊNDICE F.
102
Início: cabeceira 33 da pista do aeroporto Júlio César
Para alinhar o Aerosonde na direção ψ = 330º deve-se fornecer a atitude
inicial ao bloco do modelo não-linear que aceita um vetor na forma de quatérnions.
Para se obter a atitude em quatérnions basta utilizar as equações de (2.5) do
Capítulo 2:
q0 = ± (cos φ/2 cos θ/2 cosψ/2 + sen φ/2 sen θ/2 senψ/2)
q1 = ±(sen φ/2 cos θ/2 cosψ/2 − cos φ/2 sen θ/2 senψ/2)
q2 = ± (cos φ/2 sen θ/2 cosψ/2 + sen φ/2 cos θ/2 senψ/2)
q3 = ± (cos φ/2 cos θ/2 senψ/2 − sen φ/2 sen θ/2 cosψ/2)
com φ = 0º e θ = 0º . O quatérnion obtido foi:
[ q0
q1
q2
q3 ] = [ −0,9659 0 0 0, 2588]
T
T
Latitude: 1º 25'07, 70"S . Como o sistema de navegação trabalha com o
sistema de coordenadas em radianos utiliza-se a função dms2rad do Matlab para a
conversão, conforme:
dms2rad(-1,25,07.70);
Latitude ( rad ): -0, 02476282839003 rad
Percebe-se que o sinal de menos (-) foi colocado antes do 1 devido a coordenada
estar ao sul da linha do equador.
Longitude: 48º 27 '22, 47"W (Oeste) = -0,84572096022540 rad ;
Altitude: 1000 m ; Va = 23 m / s ; φMAX = 15º ; θ MAX = 7º .
WP 1: cabeceira 15 da pista do aeroporto Júlio César
Latitude: 1º 24 '34, 23" S = -0, 02460056125096 rad ;
Longitude: 48º 27 '55, 20"W = -0.84587963974323 rad ;
Altitude: 975 m ; Va = 23 m / s ; φMAX = 15º ; θ MAX = 7º .
WP 2: Baía do Guajará
Latitude: 1º 26 '06,59" S = -0.02504833516684 rad ;
Longitude: 48º 29 '50, 29"W = -0,84643761180881 rad ;
Altitude: 950 m ; Va = 23 m / s ; φMAX = 15º ; θ MAX = 7º .
103
WP 3: sobre a Casa das 11 Janelas, Belém-PA, Brasil
Latitude: 1º 27 '23,86" S = -0.02542295069823 rad ;
Longitude: 48º 30 '19, 66"W = -0.84658000158696 rad ;
Altitude: 1000 m ; Va = 25 m / s ; φMAX = 20º ; θ MAX = 7º .
WP 4: ponto de aproximação para a cabeceira 33
Latitude: 1º 25'31, 75" S = -0.02487942608034 rad ;
Longitude: 48º 26 '59, 06"W = -0.84560746534265 rad ;
Altitude: 1000 m ; Va = 23 m / s ; φMAX = 30º ; θ MAX = 7º .
WP 5: cabeceira 33 da pista do Aeroporto Júlio César
Latitude: 1º 25'07, 70" S = -0, 02476282839003 rad ;
Longitude: 48º 27 '22, 47"W = -0,84572096022540 rad ;
Altitude: 1000 m ; Va = 23 m / s ; φMAX = 15º ; θ MAX = 7º .
Figura 5. 15: Rota gerada pelo waypoint_gen.m.
104
5.4.1.2 Teste dos controladores com o modelo não-linear
Controlador LQG/LTR Lateral e Longitudinal com o modelo não-linear
Teste feito iniciando a simulação sob as condições do modelo linear do
Capítulo 3, solicitando ao controlador longitudinal que rastreie Varef = 23 m / s ,
θ ref = 7º e solicitando ao controlador lateral que mantenha as asas niveladas. As
respostas de θ e φ são apresentadas na Figura 5.16 e de Va em 5.17. Perceba que
θ se estabiliza em torno de 7º quando a velocidade se aproxima de 23 m / s .
Figura 5. 16: Resposta temporal do controle de
θ
e
φ
com zoom nos primeiros 4 s de simulação.
Figura 5. 17: Resposta temporal de Va para o teste dos controladores com asas niveladas.
105
Em uma segunda simulação, o controlador longitudinal permanece com as
mesmas referências, θ ref = 7º e Varef = 23 m / s , enquanto que um doublet de
φref = ±15º é aplicado na referência do controlador lateral de 0 a 30 s e de 30 a
60 s , conforme mostrado na Figura 5.18.
Figura 5. 18: Resposta temporal da simulação com doublet aplicado em
φref
e o sinal de controle.
106
O sinal de controle mostrado na Figura 5.18 representa o momento de
transição de φref quando este passa de +15º para −15º . Percebe-se que os
controladores LQG/LTR longitudinais e laterais, tal como ocorreu com a simulação
do modelo linear, envia comandos de grandes amplitudes durante um curto período
de tempo. Este comportamento já era esperado, pois estes controladores devem
comportar-se do mesmo modo nas duas simulações (com modelo linear e não
linear), a diferença é que no modelo não linear, a aeronave não responde da mesma
forma, tal como pode ser visto na Figura 5.17 a ausência do sobre-sinal de 60% na
resposta do ângulo de pitch. Outra observação, esta em relação a Figura 5.18, é
referente ao sinal de controle do acelerador, que está limitado entre 0 e 1, sendo que
no trecho visível da figura o sinal de controle do acelerador se mantém em 1 para
garantir que a velocidade de 23 m/s possa ser rastreada enquanto mantém-se o
ângulo de pitch rastreando 7º.
Com estes resultados, verificou-se que os controladores de baixo nível
projetados pelo método LQG/LTR possibilitam o controle das variáveis desejadas
com precisão estática e até com transitórios mais suaves do que havia sido
verificado no Capítulo 4. Os sinais de controle estão dentro de uma faixa aceitável e
com isso, pode-se passar ao estudo dos demais controladores das camadas mais
superiores, como os de controle de altitude e direção.
Controle fuzzy de Altitude e Direção com o modelo não-linear
Neste teste os controladores fuzzy fornecem as referências aos subcontroladores lateral e longitudinal para rastrear altitude e direção. Os limitadores
angulares estarão em φMAX = 15º e θ MAX = 7º . A simulação inicia na altitude ( h ) de
980 m e ψ = 0º para rastrear 1000 m e 45º com β → 0 e Varef = 23 m / s . Lembrando
que href = 1000 − offset , mencionado na seção 5.1 com offset = 3 . A resposta da
altitude é mostrada na Figura 5.19, φ , θ , β , e α em 5.20 e, ψ e Va em 5.21.
Percebe-se que a transição na altitude tem um valor praticamente constante nos 10
primeiros segundos, devido a tentativa de aumentar Va pela redução em θ .
107
Figura 5. 19: Resposta temporal de h e intervalo
Figura 5. 20: Resposta temporal de φ,
marcando o momento em que
θ, β, e α
θ
reduz e Va aumenta.
no controle de altitude e direção.
108
Figura 5. 21: Resposta temporal de Va e
ψ
no controle de altitude e direção.
Todas as simulações desta seção e das seções seguintes utilizam o algoritmo
de simulação ODE4 (Runge-Kutta) do Simulink com período de amostragem de
0, 02 s .
Após estes testes, verificou-se que os controladores estão garantindo bons
resultados com o modelo não-linear do Aerosonde, portanto, dá-se continuidade a
simulação da rota completa da seção 5.4.1.1.
5.4.1.3 Simulação da rota: Sobrevôo em Belém
Na Figura 5.22 é apresentado o gráfico tridimensional da rota percorrida pelo
Aerosonde em função da longitude, latitude e altitude. Na Figura 5.23 a mesma rota
é apresentada em duas dimensões em função apenas da longitude e latitude. Por
fim, a variação da altitude no tempo é mostrada na Figura 5.24, onde se percebe um
pequeno erro em regime permanente, porém, aceitável por ser menor que 60 m .
As manobras feitas pelo Aerosonde em termos de sua direção e as variações
na velocidade são apresentadas nas figuras 5.25 e 5.26, respectivamente.
109
Figura 5. 22: Gráfico 3-D da rota percorrida pelo Aerosonde UAV no Sobrevôo em Belém.
Figura 5. 23: Gráfico 2-D da rota percorrida pelo Aerosonde UAV no Sobrevôo em Belém.
110
Figura 5. 24: Variação da altitude ao longo da simulação.
Figura 5. 25: Variação de ψ ao longo da simulação.
Figura 5. 26: Variação de Va ao longo da simulação.
111
5.4.2 Simulação 2: com variação de parâmetros, perturbações e fora da região
de estabilidade programada
Nesta simulação, a rota de Sobrevôo em Belém será a mesma, exceto pela
altitude, que será 700 m mais baixa que a simulação anterior. Além desta
modificação, a massa de combustível da aeronave será reduzida em 1 kg. Estas
mudanças afetarão diretamente no comportamento aerodinâmico e na distribuição
de massa do Aerosonde.
Até o presente momento, todas as simulações estavam sendo feitas sem a
presença de ventos (constantes ou em rajadas) e turbulências, que para esta
simulação foram incluídas e possuem representação gráfica apresentada na Figura
5.27. O sistema de coordenadas adotado utiliza a convenção NED, logo, o vetor em
(5.9) apresenta ventos provenientes do sul e do oeste (devido ao sinal negativo) e
turbulências no sentido positivo (para baixo) do eixo z. O vetor em (5.9) é
multiplicado por uma sequência aleatória de média 3 , variância nula e baixa
frequência para simular as perturbações atmosféricas.
[-1 - 0,9 0, 2]T
(5.9)
Figura 5. 27: Velocidades dos ventos sul e oeste e da turbulência usada na simulação.
112
Na Figura 5.28 é apresentado o gráfico tridimensional da rota percorrida pelo
Aerosonde e, na Figura 5.29, a representação em duas dimensões da rotas com e
sem perturbações e variações de parâmetros para que seja possível comparar a
simulação 1 e 2.
A Figura 5.30 apresenta as respostas da variação da altitude e direção em
função do tempo, onde é facilmente notada a presença da turbulência afetando o
controle da altitude assim como os ventos influenciando nas correções de direção.
Apesar do erro em regime permanente na resposta da altitude, esta ainda é bem
inferior à tolerância adotada na indústria de ±60 m (Stevens e Lewis; 2003).
A velocidade sofreu diversos picos, mas manteve a média em torno de
23 m / s como pode ser visto na Figura 5.31.
Deve-se fazer uma observação em relação às perturbações atmosféricas que
foram simuladas, pois estas submeteram o Aerosonde a variações extremas,
quando na verdade, no mundo real, variações nas velocidades dos ventos não
ocorrem de forma tão abrupta como foi simulado.
Figura 5. 28: Rota 3-D percorrida pelo Aerosonde submetido a perturbações atmosféricas e variação de
parâmetros. A altitude dos waypoints entre 250 e 300 m .
113
Figura 5. 29: Gráfico 2-D da Simulação 1 e 2 juntas para comparação entre simulação com e sem perturbações e
variações de parâmetros.
Figura 5. 30: Direção e Altitude foram afetadas pelos ventos e turbulências.
114
Figura 5. 31: Va e Vamédia ao longo da simulação com perturbações e variação de parâmetros.
Com estes resultados, fica evidente que a utilização das técnicas de controle
robusto multivariável em união com o controle fuzzy, quando respeitadas as
margens de robustez e estabilidade, garantem controlar eficientemente a planta
mesmo para grandes variações no modelo e na presença de perturbações.
115
Capítulo 6
Conclusão Geral e Sugestões para novos trabalhos
No Capítulo 2, uma revisão teórica de conceitos sobre cinemática e dinâmica
de aeronaves foi apresentada, ressaltando-se as equações gerais de movimento de
corpos rígidos com seis graus de liberdade, representações de atitudes em espaços
tridimensionais por ângulos de Euler e quatérnions, e como relacionar estas
representações. Um resumo sobre aerodinâmica básica, frisando o funcionamento
de aerofólios e as principais forças e momentos envolvidos, foi também abordado.
Este estudo forneceu o conhecimento necessário para melhor compreender a
relação entre as variáveis envolvidas no problema e como são feitas simulações com
aviões.
Ainda no Capítulo 2, um estudo sobre geodésia e navegação terrestre foi
desenvolvido para fornecer subsídios ao entendimento da navegação de aeronaves
ao redor do globo terrestre e como são tratadas as informações de posicionamento
na Terra. Os modelos mais comuns do geoid foram mostrados, com noções sobre
frames de referência e sistemas de coordenadas normalmente utilizados no estudo
de veículos aéreos; mostrou-se também os eixos e ângulos utilizados, as superfícies
móveis de controle aerodinâmico dos aviões e respectivas convenções de sinais
adotadas pela indústria aeronáutica, consideradas informações de suma importância
para entender os resultados obtidos durante as simulações.
No Capítulo 3, uma apresentação detalhada dos principais blocos e funções
do Aerosim 1.2 foi fornecida, assim como o modelo do UAV Aerosonde e uma
introdução sobre veículos aéreos não tripulados; um modelo matemático linear do
Aerosonde foi obtido e analisado detalhadamente, em relação aos principais modos
naturais envolvidos na dinâmica de aviões. Esta fase de análise permitiu maior
familiarização com a ferramenta de simulação, o Aerosim, e também antecipou o
conhecimento dos principais problemas e limitações que seriam encontradas para
ajudar na solução dos problemas na fase de projeto dos controladores. Como
exemplo,
verificou-se
que
a
planta
possuía
pólos
instáveis,
apresentava
cancelamento de pólos e zeros de transmissão e relacionou-se movimentos e
comportamentos característicos da aeronave aos pólos do sistema, dando uma
interpretação física para esta fase de análise.
116
No Capítulo 4, fez-se uma revisão geral sobre os principais sistemas de
controle de vôo e de aumento de estabilidade utilizados em aviões, com
apresentação dos pilotos automáticos mais comuns e formas de avaliar o
desempenho dos controladores baseando-se em especificações de qualidade
utilizadas pela indústria aeronáutica. Fez-se também uma revisão teórica de técnicas
de controle robusto multivariável e análise no domínio da frequência, dando ênfase
sobre o regulador LQG/LTR e a aplicação deste no controle de aeronaves,
garantindo-se robustez e estabilidade em uma grande faixa de operação, ou fases
de vôo, baseadas em especificações de desempenho no domínio da freqüência.
Foram projetados os controladores lateral e longitudinal, seguindo-se à risca as
margens de robustez e estabilidade prescritas. Ainda na fase de projeto, algumas
limitações foram detectadas, estas em relação às restrições de desempenho
baseadas na análise dos ganhos principais do sistema no domínio da frequência.
Como estas restrições eram muito exigentes, tornou-se difícil a obtenção de
respostas temporais sem grandes sobre-sinais, isto para que a faixa de operação
dos controladores não fosse reduzida.
Apesar das limitações encontradas, quando os controladores LQG/LTR do
Capítulo 4 foram testados com o modelo não linear, verificou-se uma melhoria em
relação a resposta transitória, uma melhoria associada ao fato de o modelo não
linear considerar a dinâmica dos atuadores do Aerosonde.
No Capítulo 5, foi finalizado o sistema de controle de vôo do Aerosonde com o
projeto de controladores fuzzy, do tipo Mamdani, para altitude e direção, baseado
em um resumo teórico sobre o tipo de controlador utilizado; apresentou-se o sistema
de navegação projetado e como este determina a trajetória baseando-se nas
informações dos waypoints. O piloto automático completo do Aerosonde foi
apresentado bloco por bloco assim como a estrutura interna utilizada; simulações
individuais com o modelo não linear do Aerosonde, com 6-DOF, foram feitas com
cada bloco desenvolvido, culminando em testes finais com planejamento de missão,
criação dos waypoints e simulação do sistema em duas situações, sem e com
perturbações e variação de parâmetros. Nas duas simulações o Aerosonde
conseguiu percorrer a rota desejada tendo bons resultados mesmo quando testado
fora do ponto de operação ao qual os controladores foram projetados, o que indica
que a fusão das técnicas inteligentes e modernas neste tipo de aplicação pode servir
117
como alternativa para aeronaves que não exijam uma grande faixa de operação em
termos de altitude e velocidade.
Acredita-se que os resultados aqui apresentados sejam de motivação para
novos estudos de técnicas de controle aplicadas a veículos aéreos, pois o método
adotado para controlar a planta em um envelope altitude-velocidade razoavelmente
grande baseou-se nos projetos realizados com um único modelo linear submetido a
uma forte restrição para a garantia de robustez e estabilidade. Logo, tentar explorar
outras situações mais amplas, com múltiplos modelos locais e utilização de gainscheduling, poderia manter ou ampliar o envelope de operação, com maior liberdade
ao projetista para moldar as respostas temporais. Ainda no âmbito de se trabalhar
com múltiplos modelos lineares locais, seria interessante experimentar outras
técnicas de fusão de controle inteligente e controle moderno, utilizando, por
exemplo, um Compensador Paralelo Distribuído e verificação da estabilidade do
modelo global pelo segundo método de Lyapunov, com solução de Desigualdades
Matriciais Lineares (Passino e Yurkovich; 1998).
Para uma alternativa mais simples, sugere-se manter a mesma estrutura de
controladores abordada neste trabalho, mas com alterações nos controladores fuzzy,
incluindo novas variáveis de entrada, como por exemplo, a derivada do erro, e
aumentar o número de conjuntos fuzzy. Outras linhas de estudo poderiam ser
seguidas também no ramo dos sistemas de navegação, fornecendo maior
funcionalidade ao sistema de obtenção de trajetória, como a inclusão da altitude nas
restrições de waypoints e obtenção de uma trajetória de subida (rampa de subida)
pelo método de triangulação, com a possibilidade de determinação de uma razão de
subida ótima baseando-se na economia de combustível.
Por fim, seria interessante testar o piloto automático em um protótipo real,
como em um aeromodelo adequado.
118
Referências Bibliográficas
ANDERSON, B.D.O., MOORE, J.B. . Optimal Control: Linear Quadratic Methods.
Prentice-Hall International, Inc. . 1989.
CAMPOS, M., ISCOLD, P., TORRES, L.A.B., AGUIRRE, L.A. . SiDeVAAN:
Simulação e Desenvolvimento de Veículos Aéreos Autônomos e NãoTripulados. VIII Simpósio Brasileiro de Automação Inteligente – SBAI. 2007.
Florianópolis-SC.
CAVALCANTE, C.M.C. . Sistema de Navegação para Helicópteros Não Tripulados
Utilizando Controlador Nebuloso. Dissertação de Mestrado submetida ao curso
de Engenharia Elétrica da Universidade Federal de Santa Catarina – UFSC.
Florianópolis-SC, 1994.
DOITSIDIS, L., VALAVANIS, K.P., TSOURVELOUDIS, N.C., KONTITSIS, M. . A
Framework for Fuzzy Logic Based UAV. Proceedings of the 2004 IEEE
International Conference on Robotics & Automation. New Orleans, LA. April
2004.
DOYLE, J.C., STEIN, G. . Robustness with Observers. IEEE Transactions on
Automatic Control, vol. AC-24, no. 4, Agosto, 1979, pp. 607-611.
ETKIN, B. . Dynamics of Flight – Stability and Control. John Wiley & Sons. 1959.
GREEN, M., LIMEBEER, D.J.N. . Linear Robust Control. Prentice Hall, New Jersey,
U.S.A. 1995.
HOOVER, R.C. . Fusion of Hard and Soft Computing for Guidance, Navigation and
Control of Uninhabited Aerial Vehicles. Thesis submitted in partial fulfillment of
the Requirements for the degree of M.Sc. in Measurement and Control
Engineering. Idaho State University, Dec. 2004.
Introduction
to
UAV
Systems
–
UAV
Center
Co.
Ltd..
Disponível
em:
<http://www.uavcenter.com>. Acesso em: 30 mar. 2007.
Jeppesen General Navigation Manual – Jeppesen Sanderson, Inc.. 2005. Disponível
em: <http://www.jeppesen.com>. Acesso em: 30 mar. 2007.
LY, Uy-Loi. Stability and Control of Flight Vehicle. Department of Aeronautics and
Astronautics – University of Washington. September 29, 1997.
119
MILEVA, B., BOSKOSKI, P., DESKOSKI, S. . Fuzzy Logic Controlo f the Hight of the
Airplane. 6th International PhD Workshop on Systems and Control, October 4-8,
2005 Izola, Slovenia.
NERIS, L.O. . Um piloto automático para as aeronaves do projeto ARARA.
Dissertação de mestrado defendida no Instituto de Ciências Matemáticas e de
Computação, USP, 2001.
PASSINO, K.M., YURKOVICH, S. . Fuzzy Control. Addison Wesley Longman, Inc.
1998.
Precision
Manuals
Boeing
737NG
documentation.
Disponível
em:
<http://www.precisionmanuals.com>. Acesso em: 30 mar. 2007.
RAMOS, J., D’ABREU, J., NEVES, O., FIGUEIREDO, D., TANURE, L., HOLANDA,
F. . Iniciativa para Robótica Pedagógica Aberta e de Baixo Custo para Inclusão
Social e Digital no Brasil. VIII Simpósio Brasileiro de Automação Inteligente –
SBAI. 2007. Florianópolis-SC.
ROSKAM, J. . Airplane Flight Dynamics and Automatic Flight Controls. University of
Kansas, Lawrence, Kansas. Copyright of Roskam Aviation and Engineering
Corporation. 1979.
SKOGESTAD, S., POSTLETHWAITE, I. . Multivariable Feedback Control: Analysis
and design. John Wiley & Sons. 1996.
SILVEIRA, A.S., LOPES, A.M.O., da COSTA Jr, C.T., BARREIROS, J.A.L. . Sistema
de Controle fuzzy para Veículo Aéreo Não Tripulado. VIII Simpósio Brasileiro de
Automação Inteligente – SBAI. 2007. Florianópolis-SC.
STEVENS, B.L. e LEWIS, F.L. . Aircraft Control and Simulation. 2ª Ed. 2003. John
Wiley & Sons.
SUGENO, M., HIRANO, I., NAKAMURA, S., KOTSU, S. . Development of na
Intelligent Unmanned Helicopter. Tokyo Institute of Technology, Japan. IEEE
1995.
WANG, Li-Chin. A course in fuzzy systems and control. Prentice Hall PTR. 1ª Ed.
1997.
ZAREI, J., MONTAZERI, A., MOTLAGH, M.R.J., POSHTAN, J. . Design and
comparison of LQG/LTR and H ∞ controllers for a VSTOL flight control system.
Journal of the Franklin Institute. 2006.
120
APÊNDICE A - Instalação do Aerosim 1.2 para MATLAB 6 ou superior
Para instalar o Aerosim Blockset no Matlab, execute o arquivo setup.exe que
acompanha o pacote de instalação que pode ser encontrado em www.udynamics.com. Por padrão, após a instalação, os arquivos do Aerosim ficarão na
pasta c:\Arquivos de Programas\Aerosim ou o usuário pode selecionar outro diretório
que deverá ser conhecido para uso em outro procedimento.
Após a instalação, o arquivo c:\Pasta_do_Matlab\toolbox\local\pathdef.m será
modificado automaticamente para a inclusão da pasta do Aerosim ao Matlab. Para
verificar se a instalação foi bem sucedida, abra o Matlab e o Simulink. Dentro do
Simulink, procure na lista de blocos pelo item AeroSim Blockset. Se não aparecer
esta opção, será necessário atualizar o path do Matlab manualmente. Feche o
Simulink e, na janela principal do Matlab, vá ao Menu arquivo e clique sobre Set
Path; dê outro clique sobre o botão “Add with subfolders” e adicione a pasta onde
o Aerosim 1.2 foi instalado. Feche e abra o Matlab e vá novamente ao Simulink para
verificar se o AeroSim Blockset aparece na lista de itens.
Conectando o Aerosim ao Microsoft Flight Simulator
Primeiro deve-se obter o software FSUIPC de Peter Dowson na Internet em
www.schiratti.com. Para instalá-lo no Flight Simulator, copie o arquivo fsuipc.dll para
a pasta c:\Seu_Diretório_do_FlightSim\Modules.
No Matlab, prepare seu diagrama de simulação no Simulink incluindo o bloco
\AeroSim Blockset\Pilot Interface\FS Interface. Este bloco deverá receber a
informação dos sensores do bloco da planta, contendo: posição (latitude, longitude e
altitude), atitude em ângulos de Euler (radianos) e a velocidade do fluxo de ar.
Dando um duplo clique sobre o bloco FS Interface, deve-se especificar o período de
amostragem, devendo este ser menor ou igual ao que se está utilizando na
simulação do Simulink. Um exemplo pronto encontra-se no pacote de arquivos do
Aerosim. Para executá-lo, basta digitar aerosonde_demo_msfs na janela de
comandos do Matlab.
121
Dependendo da demanda de processamento para efetuar a simulação no
Matlab/Simulink e apresentar o andamento em tempo real no Flight Simulator, pode
ser necessário utilizar dois computadores. Caso este seja o caso, o usuário
precisará do software Wide FS encontrado também na página de Peter Dowson. O
Wide FS é dividido em dois softwares, um cliente e um servidor, que operam em
uma rede TCP/IP. Para utilizá-lo, basta copiar o arquivo servidor no mesmo
computador onde o Flight Simulator será executado e, o cliente, onde o
Matlab/Simulink estará rodando. No Wide FS Cliente, deve-se informar o endereço
IP do servidor para que os dados da simulação sejam passados via rede ao
computador com o Flight Simulator.
Para maiores informações, recomenda-se visitar a página de Peter Dowson e
do Aerosim na Internet.
122
APÊNDICE B – Controlador LQG/LTR Lateral do Aerosonde
close all;
%clear;
% Lateral-directional Dynamics
% -----------------------------% State vector: x = [v p r phi psi]
% Input vector: u = [aileron rudder]
% Output vector: y = [beta p r phi psi]
% State matrix: A =
A=[-0.5895
1.7312 -22.9345
9.7991
0;
9.1681
0
0;
0.6278 -2.4709 -0.9582
0
0;
0.0755
-0.0000
0;
1.0028
-0.0000
0];
-3.8720 -19.0490
0
1.0000
0
0
% Control matrix: B =
B=[-1.1552
2.9486;
-101.4284
1.8250;
-3.9992 -18.6309;
0
0;
0
0];
% Observation matrix: C =
C=[ 0
0
0.0435 0
0
1
0
0
0;
0];
% Plotando SVs do sistema original
original_sys=ss(A,B,C,0);
w=logspace(-3,3,100);
figure(1);
sigma(original_sys,w), grid, title('SVs do sistema original'), hold on
% Polos, fator de amortecimento e frequencias
disp('Eigenvalues & dampings of the original system');
damp(original_sys);
% Aumentando o sistema original com pre-comp. e integradores nas entradas
% Matriz de ganho DC inversa
P=inv(dcgain(original_sys));
123
Aint=-0.001*diag([1 1]); % integradores em -0.001 prox. a origem
% indice 'a' para sistema augmented ou aumentado
Aa=[A B; zeros(2,5) Aint];
Ba=[0*B; P];
Ca=[C zeros(2)];
Da=[zeros(2,2)];
augmented_sys=ss(Aa,Ba,Ca,Da);
figure(1);
sigma(augmented_sys,w), title('SVs do sistema aumentado com P'), legend('original
SVs','augmented SVs')
% Polos, fator de amortecimento e frequencias do sistema aumentado
disp('Eigenvalues & dampings of the augmented system');
damp(augmented_sys);
% Incertezas de alta frequencia
high_bound=tf([20],[1 2]);
sigma(high_bound,w);
% Filtro de Kalman
%
[ v
p
r phi psi
i1 i2]
Qk=diag([1 1 1 1 1 0.1 0.1]); % Matriz de pesos Q do filtro
rf=1e+2; % rho
Rk=eye(2); % Matriz de pesos R do filtro
disp('Kalman gain L');
L=lqr(Aa',Ca',rf*Qk,Rk)' % Ganho otimo de Kalman
kalman_sys=ss(Aa,L,Ca,Da); % Eq. de estados do filtro
figure(2);
sigma(kalman_sys,w), title('SVs com Filtro de Kalman'), grid;
% Margem de ganho e fase
alfa0=sigma(kalman_sys,w,2);
alfa0=min(alfa0(1,:));
MG=1/(1-alfa0);
MF=2*asin(alfa0/2)*180/pi;
% Testa alcancabilidade & observabilidade
Alc = obsv(Aa,sqrt(Qk));
unAlc = length(Aa)-rank(Alc);
124
if unAlc == 0
disp('alcancavel');
else
disp('nao alcancavel');
end
Ob=obsv(Aa',Ba');
unob=length(Aa)-rank(Ob);
if unob == 0
disp('observavel');
else
disp('nao observavel');
end
% Resposta ao degrau com Filtro de Kalman
t=[0:0.01:5];
y=step(Aa-L*Ca,L,Ca,Da,1,t); % mudar de 1 para 2 para testar degrau na 2a entrada
figure(3) % CHANNEL 1
plot(t,y), grid, title('Resposta ao degrau aplicado na referencia 1'), legend(' ',' ');
% polos, fator amort. e frequencias com (Aa-LCa)
disp('Eigenvalues & dampings of the Kalman system (Aa-LCa)');
damp(ss(Aa-L*Ca,L,Ca,Da));
% Regulador LQR
%
[ v
p
r phi psi
i1 i2]
rc=1.000e-2; %rho
Qlq=diag([1e2 1e0 1e0 1e6 1e-6
1 1]); % Matriz de pesos Q do LQR
Rlq=rc*eye(2); % Matriz de pesos R do LQR
disp('LQR gain Klq');
Klq=lqr(Aa,Ba,Qlq,Rlq) % Ganho otimo do LQR
lqr_sys=ss(Aa,Ba,Klq,Da); % Eq. estados do LQR
figure(4);
sigma(lqr_sys,w), title('SVs com LQR'), grid;
% LQG/LTR Regulator
Alqg=[Aa -Ba*Klq; zeros(7) Aa-L*Ca-Ba*Klq]; % Sistema em MA com LQG/LTR
Blqg=[zeros(7,2); -L];
Clqg=[Ca zeros(2,7)];
Dlqg=[zeros(2)];
125
lqg_sys=ss(Alqg,Blqg,Clqg,Dlqg);
figure(4);
sigma(lqg_sys,w), title('Valores Singulares com LQG/LTR'), grid;
% Margem de ganho e fase
alfa0=sigma(lqg_sys,w,2);
alfa0=min(alfa0(1,:));
MG=1/(1-alfa0);
MF=2*asin(alfa0/2)*180/pi;
% polos, fator amort. e freq. em MF com LQG/LTR
disp('Eigenvalues & dampings of the system with LQG/LTR');
damp(ss(Alqg-Blqg*Clqg,Blqg,Clqg,Dlqg));
% Resposta ao degrau em MF com LQG/LTR
t=[0:0.01:10];
y=step(Alqg-Blqg*Clqg,Blqg,Clqg,Dlqg,1,t); % mudar para 2 para aplicar degrau na 2a
entrada
figure(5) % CHANNEL 1
plot(t,y), grid, title('Resposta ao degrau aplicado na referencia phi com LQG/LTR'),
legend(' ',' ');
pause
close all
126
APÊNDICE C – Controlador LQG/LTR Longitudinal do Aerosonde
close all;
%clear;
% Longitudinal Dynamics
% ----------------------% State vector: x = [u w q theta h Omega]
% Input vector: u = [elevator throttle]
% Output vector: y = [Va alpha q theta h]
% State matrix: LongA =
LongA=[-0.2029
0.6110 -1.7044 -9.7991 -0.0001
0.0100;
-0.5670 -3.8086 22.4291 -0.7397
0.0010
0.4906 -4.2213 -4.3901
0
0.0000 -0.0078;
0
0
0
1.0000
0.0753
-0.9972
0
29.7630
2.2466
0
22.9997
0
0;
0
0;
0
0;
-0.0379 -2.6886];
% Control matrix: LongB =
LongB=[ 0.3132
0;
-1.9847
0;
-27.5486
0;
0
0;
0
0;
0 349.2659];
% Observation matrix: LongC =
LongC=[ 0.9972
0
0.0753
0
0
0
0;
0
0
1
0
0];
% Plotando SVs do sistema original
original_sys=ss(LongA,LongB,LongC,0);
w=logspace(-3,3,100);
figure(1);
sigma(original_sys,w), grid, title('SVs do sistema original'), hold on
% polos, fator de amortecimento e frequencias do sistema original
disp('Eigenvalues & dampings of the original system');
damp(original_sys);
127
% Aumentando o sistema com pre-comp. e integradores na entrada
% Matriz do ganho DC inversa
LongP=inv(dcgain(original_sys));
LongAint=-0.001*diag([1 1]); % integradores em -0.001
% indice 'a' para sistema aumentado ou augmented
LongAa=[LongA LongB; zeros(2,6) LongAint];
LongBa=[0*LongB; LongP];
LongCa=[LongC zeros(2)];
LongDa=[zeros(2,2)];
augmented_sys=ss(LongAa,LongBa,LongCa,LongDa); % sistema aumentado
figure(1);
sigma(augmented_sys,w), title('SVs do sistema aumentado com pre-compensador'),
legend('original SVs','augmented SVs')
% polos, fator de amort., frequencias do sistema aumentado
disp('Eigenvalues & dampings of the augmented system');
damp(augmented_sys);
% Incertezas de alta frequencia:
high_bound=tf([20],[1 2]);
sigma(high_bound,w);
% Filtro de Kalman
%
[u w q theta h Omega int1 int2]
LongQk=diag([1 1 1 1 1 1 30 90]); % Matriz de pesos Q do filtro
Longrf=1e1; %rho
LongRk=eye(2); % Matriz de pesos R do filtro
disp('Kalman gain LongL');
LongL=lqr(LongAa',LongCa',Longrf*LongQk,LongRk)' % Ganho otimo de Kalman
kalman_sys=ss(LongAa,LongL,LongCa,LongDa);
figure(2);
sigma(kalman_sys,w), title('SVs com Filtro de Kalman'), grid;
% Margem de ganho e fase
alfa0=sigma(kalman_sys,w,2);
alfa0=min(alfa0(1,:));
MG=1/(1-alfa0);
MF=2*asin(alfa0/2)*180/pi;
128
% Testa alcancabilidade & observabilidade
Alc = obsv(LongAa,sqrt(LongQk));
unAlc = length(LongAa)-rank(Alc);
if unAlc == 0
disp('alcancavel');
else
disp('nao alcancavel');
end
Ob=obsv(LongAa',LongBa');
unob=length(LongAa)-rank(Ob);
if unob == 0
disp('observavel');
else
disp('nao observavel');
end
% Plotando a resposta ao degrau com o Filtro de Kalman
t=[0:0.01:5];
y=step(LongAa-LongL*LongCa,LongL,LongCa,LongDa,1,t);
figure(3) % CHANNEL 1
plot(t,y), grid, title('Resposta ao degrau aplicado em Rtheta do sistema em M.F. com
Filtro de Kalman'), legend('
','
');
t=[0:0.01:5];
y=step(LongAa-LongL*LongCa,LongL,LongCa,LongDa,2,t);
figure(3) % CHANNEL 2
plot(t,y), grid, title('Resposta ao degrau aplicado em Rva do sistema em M.F. com
Filtro de Kalman'), legend('
','
');
% polos, amortecimento, frequencia com (LongAa-LongLLongCa)
disp('Eigenvalues & dampings of the Kalman system (LongAa-LongLLongCa)');
damp(ss(LongAa-LongL*LongCa,LongL,LongCa,LongDa));
% Regulador LQR
%
[u w q theta h Omega int1 int2]
Longrc=1e-2; %rho
LongQlq=diag([5e3 1e1 1e0 1e6 1e-3 1e-3
1e-3 1e-3]); % Matriz de pesos Q do LQR
LongRlq=Longrc*eye(2); % Matriz de pesos R do LQR
disp('LQR gain LongKlq');
129
LongKlq=lqr(LongAa,LongBa,LongQlq,LongRlq) % ganho otimo do LQR
lqr_sys=ss(LongAa,LongBa,LongKlq,LongDa);
figure(4);
sigma(lqr_sys,w), title('SVs with LQR'), grid;
% Regulador LQG/LTR
LongAlqg=[LongAa -LongBa*LongKlq; zeros(8) LongAa-LongL*LongCa-LongBa*LongKlq];
LongBlqg=[zeros(8,2); -LongL];
LongClqg=[LongCa zeros(2,8)];
LongDlqg=[zeros(2)];
lqg_sys=ss(LongAlqg,LongBlqg,LongClqg,LongDlqg);
figure(4);
sigma(lqg_sys,w), title('SVs com o controlador LQG/LTR'), grid;
% Margem de ganho e fase
alfa0=sigma(lqg_sys,w,2);
alfa0=min(alfa0(1,:));
MG=1/(1-alfa0);
MF=2*asin(alfa0/2)*180/pi;
% polos, amortecimento, frequencia com LQG/LTR
disp('Eigenvalues & dampings of the system with LQG/LTR');
damp(ss(LongAlqg-LongBlqg*LongClqg,LongBlqg,LongClqg,LongDlqg));
% Resposta ao degrau em MF com LQG/LTR
t=[0:0.01:10];
y=step(LongAlqg-LongBlqg*LongClqg,LongBlqg,LongClqg,LongDlqg,1,t);
figure(5) % CHANNEL 1
plot(t,y), grid, title('Resposta ao degrau aplicado em Rtheta com LQG/LTR em M.F.'),
legend('
','
');
t=[0:0.01:10];
y=step(LongAlqg-LongBlqg*LongClqg,LongBlqg,LongClqg,LongDlqg,2,t);
figure(5) % CHANNEL 2
plot(t,y), grid, title('Resposta ao degrau aplicado em Rva com LQG/LTR em M.F.'),
legend('
','
');
pause
close all
130
% Inicializa variaveis para simulacao conjunta com o Simulink
global next_wp lat_d lon_d spd_ref h_ref theta_max phi_max;
% Waypoint 1: dados devem ser substituidos de acordo com o interesse da
% missao. As variaveis abaixo foram carregadas para o primeiro waypoint da
% simulacao intitulada Sobrevoo em Belem apresentada no capitulo 5.
lat_d=-0.02460056125096;
lon_d=-0.84587963974323;
h_ref=975;
spd_ref=23;
theta_max=7;
phi_max=15;
next_wp=1;
131
APÊNDICE D – Controlador fuzzy de Altitude
function [u]=FLC_Altitude(erro);
e=erro(1); % erro de altitude
%iniciando as variaveis
Mng=0;Mz=0;Mpg=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Pertinencias para o erro
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if e <= -20
Mng=1; % pertinencia erro negativo
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>-20 & e<0
% funcao trapezoidal negativo
b=-3000;a=b;
c=-20;
d=0;
Mng=(d-e)/(d-c); % pertinencia erro negativo
end
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>-20 & e<0
%funcao triangular zero
a=-20;
c=0;b=c;
d=20;
Mz=(e-a)/(b-a); % Pertinencia zero
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e==0
Mz=1; % Pertinencia zero
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>0 & e<20
%Zero
a=-20;
c=0;b=c;
d=20;
Mz=(d-e)/(d-c); % Pertinencia zero
132
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>0 & e<20
% Positivo
a=0;
b=20;
d=3000;c=d;
Mpg=(e-a)/(b-a); % Pertinencia positivo
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>=20
Mpg=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Defuzzificacao por media dos centros
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M(1)=min(Mng); % Pertinencia negativo
B(1)=0.25; % base de DIMINUIR no univ. discurso de saida
M(2)=min(Mz); % Pertinencia zero
B(2)=0.5; % base de NEUTRALIZAR no univ. discurso de saida
M(3)=min(Mpg); % Pertinencia positivo
B(3)=0.75; % base de AUMENTAR no univ. discurso de saida
num=0;den=0;
for i=1:3,
num=num+B(i)*M(i);
den=den+M(i);
end
u=num/den; % referencia pitch saida
133
APÊNDICE E – Controlador fuzzy de Direção
function [u]=FLC_Heading(erro);
e=erro(1); % erro de direcao
if abs(e) > 180 % trata o sinal do erro
e = (abs(e)-180) - 180;
end
%iniciando as variaveis
Mng=0;Mz=0;Mpg=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Pertinencias para o erro
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if e <= -90.5
Mng=1; % Pertinencia Negativo
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>-90.5 & e<0
% Funcao trapezoidal Negativo
b=-360;a=b;
c=-90.5;
d=0;
Mng=(d-e)/(d-c);
end
%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>-90.5 & e<0
% Funcao triangular Zero
a=-90.5;
c=0;b=c;
d=90.5;
Mz=(e-a)/(b-a); % Pertinencia Zero
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e==0
Mz=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>0 & e<90.5
134
% Funcao Triangular Zero
a=-90.5;
c=0;b=c;
d=90.5;
Mz=(d-e)/(d-c);
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>0 & e<90.5
% Funcao Trapezoidal Positivo
a=0;
b=90.5;
d=360;c=d;
Mpg=(e-a)/(b-a); % Pertinencia Positivo
end
%%%%%%%%%%%%%%%%%%%%%%%%%
if e>=90.5
Mpg=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Defuzzificacao por media dos centros
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M(1)=min(Mng); % pert. negativo
B(1)=-0.5; % base funcao ESQUERDA
M(2)=min(Mz); % pert. zero
B(2)=0; % base funcao NIVELAR
M(3)=min(Mpg); % pert. positivo
B(3)=0.5; % base funcao DIREITA
num=0;den=0;
for i=1:3,
num=num+B(i)*M(i);
den=den+M(i);
end
u=num/den; % saida: referencia angulo de bank
135
APÊNDICE F – Gerador de waypoints
% Gerador de waypoints - visualiza rota da missao
clear waypoint;
waypoint(:,:,1) = [-0.02460056125096,-0.84587963974323,975];
waypoint(:,:,2) = [-0.02504833516684,-0.84643761180881,950];
waypoint(:,:,3) = [-0.02542295069823,-0.84658000158696,1000];
waypoint(:,:,4) = [-0.02487942608034,-0.84560746534265,1000];
waypoint(:,:,5) = [-0.02476282839003,-0.84572096022540,1000];
%waypoint(:,:,n) = [Lat,Lon,Alt]; Especifica mais waypoints
hold on
[tmp1,tmp2,wp_number]=size(waypoint);
for i=1:wp_number % plota os waypoints em tres dimensoes
wp=waypoint(:,:,i);
plot3(wp(2),wp(1),wp(3),'r:o','LineWidth',3),title('waypoints')...
,xlabel('Longitude (rad)'),ylabel('Latitude (rad)'),zlabel('Altitude (m)'),grid on
text(wp(2),wp(1),wp(3),[' WP ' num2str(i)])
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
% DMS2RAD Converte angulos a partir de grau:minuto:segundo para radianos
%
%
rad = dms2rad(grau,minuto,segundo)
%
% Minutos e segundos devem estar entre 0 e 60
136
APÊNDICE G – Sistema de Navegação
function [psi_d]=psi_waypoint(input); %calcula trajetoria
psi_uav=input(1); % direcao atual
beta_uav=input(2); % sideslip atual
lat_uav=input(3); % latitude atual fornecida pelo GPS
lon_uav=input(4); % longitude atual fornecida pelo GPS
% define variaveis globais para afetarem o Workspace do Matlab
global lat_d lon_d next_wp h_ref spd_ref theta_max phi_max;
if next_wp < 5 % verifica se todos os waypoints foram concluidos
% caso ainda hajam waypoints, compara (lat,lon) do waypoint
% (lat_d,lon_d) com as respectivas atuais do GPS.
if ( (abs(lat_d-lat_uav) < 0.00001) & (abs(lon_d-lon_uav) < 0.00001) )
next_wp=next_wp+1; % atualiza para o proximo waypoint
end
end
switch next_wp % switch case para setar globalmente as novas variaveis
case 2
lat_d=-0.02504833516684;
lon_d=-0.84643761180881;
h_ref=950;
spd_ref=23;
theta_max=7;
phi_max=15;
case 3
lat_d=-0.02542295069823;
lon_d=-0.84658000158696;
h_ref=1000;
spd_ref=23;
theta_max=7;
phi_max=20;
case 4
lat_d=-0.02487942608034;
lon_d=-0.84560746534265;
h_ref=1000;
spd_ref=23;
theta_max=7;
137
phi_max=30;
case 5
lat_d=-0.02476282839003;
lon_d=-0.84572096022540;
h_ref=1000;
spd_ref=23;
theta_max=7;
phi_max=15;
%
case n
% n pertence aos inteiros, para cada waypoint desejado
%
lat_d;
%
lon_d=;
%
h_ref=;
%
spd_ref=;
%
theta_max=;
%
phi_max=;
end
psi_wp=atan2((lon_d-lon_uav),(lat_d-lat_uav))*(180/pi); % calcula trajetoria
if psi_wp < 0 % trata psi para entregar os dados entre 0 e 359 (graus)
psi_wp = 360 - abs(psi_wp);
end
psi_d=[psi_wp h_ref-3 spd_ref next_wp]; % vetor de saidas da funcao
% (h_ref-3 corresponde ao offset mencionado no capitulo 5)
138
ANEXO A – Software de Configuração do Aerosonde UAV
% AIRCRAFT CONFIGURATION SCRIPT
% Aerosonde UAV - sample model from AeroSim Library
% Copyright 2002 Unmanned Dynamics, LLC
% Revision: 1.0 Date: 05/13/2002
% Clear workspace
clear all;
% Name of the MAT-file that will be generated
cfgmatfile = 'aerosondecfg';
%%% AERODYNAMICS %%%
% Aerodynamic force application point (usually the aerodynamic center)[x y z]
rAC = [0.1425 0 0]; % m
%%% Aerodynamic parameter bounds %%%
% Airspeed bounds
VaBnd = [15 50]; % m/s
% Sideslip angle bounds
BetaBnd = [-0.5 0.5]; % rad
% Angle of attack bounds
AlphaBnd = [-0.1 0.3]; % rad
%%% Aerodynamic reference parameters %%%
% Mean aerodynamic chord
MAC = 0.189941; % m
% Wind span
b = 2.8956; % m
% Wing area
S = 0.55; % m^2
% ALL aerodynamics derivatives are per radian:
%%% Lift coefficient %%%
% Zero-alpha lift
CL0 = 0.23;
% alpha derivative
CLa = 5.6106;
139
% Lift control (flap) derivative
CLdf = 0.74;
% Pitch control (elevator) derivative
CLde = 0.13;
% alpha-dot derivative
CLalphadot = 1.9724;
% Pitch rate derivative
CLq = 7.9543;
% Mach number derivative
CLM = 0;
%%% Drag coefficient %%%
% Lift at minimum drag
CLmind = 0.23;
% Minimum drag
CDmin = 0.0434;
% Lift control (flap) derivative
CDdf = 0.1467;
% Pitch control (elevator) derivative
CDde = 0.0135;
% Roll control (aileron) derivative
CDda = 0.0302;
% Yaw control (rudder) derivative
CDdr = 0.0303;
% Mach number derivative
CDM = 0;
% Oswald's coefficient
osw = 0.75;
%%% Side force coefficient %%%
% Sideslip derivative
CYbeta = -0.83;
% Roll control derivative
CYda = -0.075;
% Yaw control derivative
CYdr = 0.1914;
% Roll rate derivative
CYp = 0;
% Yaw rate derivative
CYr = 0;
140
%%% Pitch moment coefficient %%%
% Zero-alpha pitch
Cm0 = 0.135;
% alpha derivative
Cma = -2.7397;
% Lift control derivative
Cmdf = 0.0467;
% Pitch control derivative
Cmde = -0.9918;
% alpha_dot derivative
Cmalphadot = -10.3796;
% Pitch rate derivative
Cmq = -38.2067;
% Mach number derivative
CmM = 0;
%%% Roll moment coefficient %%%
% Sideslip derivative
Clbeta = -0.13;
% Roll control derivative
Clda = -0.1695;
% Yaw control derivative
Cldr = 0.0024;
% Roll rate derivative
Clp = -0.5051;
% Yaw rate derivative
Clr = 0.2519;
%%% Yaw moment coefficient %%%
% Sideslip derivative
Cnbeta = 0.0726;
% Roll control derivative
Cnda = 0.0108;
% Yaw control derivative
Cndr = -0.0693;
% Roll rate derivative
Cnp = -0.069;
% Yaw rate derivative
Cnr = -0.0946;
141
%%% PROPELLER %%%
%Propulsion force application point (usually propeller hub) [x y z]
rHub = [0 0 0]; % m
% Advance ratio vector
J = [-1 0 0.1 0.2 0.3 0.35 0.4 0.45 0.5 0.6 0.7 0.8 0.9 1 1.2 2];
% Coefficient of thrust look-up table CT = CT(J)
CT = [0.0492 0.0286 0.0266 0.0232 0.0343 0.034 0.0372 0.0314 0.0254 0.0117 -0.005 -0.0156 0.0203 -0.0295 -0.04 -0.1115];
% Coefficient of power look-up table CP = CP(J)
CP = [0.0199 0.0207 0.0191 0.0169 0.0217 0.0223 0.0254 0.0235 0.0212 0.0146 0.0038 -0.005 0.0097 -0.018 -0.0273 -0.0737];
% Propeller radius
Rprop = 0.254; % m
% Propeller moment of inertia
Jprop = 0.002; % kg*m^2
%%% ENGINE %%%
% Engine rpm vector
RPM = [1500 2100 2800 3500 4500 5100 5500 6000 7000]; % rot per min
% Manifold pressure vector
MAP = [60 70 80 90 92 94 96 98 100]; % kPa
% Sea-level fuel flow look-up table fflow = fflow(RPM, MAP)
% RPM -> rows, MAP -> columns
FuelFlow = [
31 32 46 53 55 57 65 73 82
40 44 54 69 74 80 92 103 111
50 63 69 92 95 98 126 145 153
66 75 87 110 117 127 150 175 190
83 98 115 143 148 162 191 232 246
93 102 130 159 167 182 208 260 310
100 118 137 169 178 190 232 287 313
104 126 151 184 191 206 253 326 337
123 144 174 210 217 244 321 400 408
]; % g/hr
% Sea-level power look-up table P = P(RPM, MAP)
% RPM -> rows, MAP -> columns
142
Power = [
18.85 47.12 65.97 67.54 69.12 67.54 67.54 69.12 86.39
59.38 98.96 127.55 149.54 151.74 160.54 178.13 200.12 224.31
93.83 149.54 187.66 237.5 249.23 255.1 307.88 366.52 398.77
109.96 161.27 245.57 307.88 326.2 351.86 421.5 491.14 531.45
164.93 245.04 339.29 438.25 447.68 494.8 565.49 673.87 772.83
181.58 245.67 389.87 496.69 528.73 571.46 662.25 822.47 993.37
184.31 293.74 403.17 535.64 570.2 622.04 748.75 956.09 1059.76
163.36 276.46 420.97 565.49 609.47 691.15 860.8 1130.97 1193.81
124.62 249.23 417.83 586.43 645.07 762.36 996.93 1246.17 1429.42
]; % W
% Sea-level pressure and temperature at which the data above is given
pSL = 102300; % Pa
TSL = 291.15; % deg K
% Engine shaft moment of inertia
Jeng = 0.0001; % kg*m^2
%%% INERTIA %%%
% Empty aircraft mass (zero-fuel)
mempty = 8.5; % kg
% Gross aircraft mass (full fuel tank)
mgross = 13.5; % kg
% Empty CG location [x y z]
CGempty = [0.156 0 0.079]; % m
% Gross CG location [x y z]
CGgross = [0.159 0 0.090]; % m
% Empty moments of inertia [Jx Jy Jz Jxz]
Jempty = [0.7795 1.122 1.752 0.1211]; % kg*m^2
% Gross moments of inertia [Jx Jy Jz Jxz]
Jgross = [0.8244 1.135 1.759 0.1204]; % kg*m^2
%%% OTHER SIMULATION PARAMETERS %%%
% WMM-2000 date [day month year]
dmy = [13 05 2002];
% Save workspace variables to MAT file
save(cfgmatfile);
143
% Output a message to the screen
fprintf(strcat('\n Aircraft configuration saved as:\t', strcat(cfgmatfile),'.mat'));
fprintf('\n');
144
ANEXO B – Software de Trim e Linearização do Aerosonde
% trim_aerosonde
% A Matlab script that trims the nonlinear aircraft model for a chosen
% flight condition and extracts the aircraft linear model
%
% Unmanned Dynamics, LLC
% October 1, 2002
%
% The trim parameters structure has the following components:
% 1. Simulation settings
% TrimParam.SampleTime = simulation sample time
% TrimParam.FinalTime = simulation final time (used only for trim function)
% TrimParam.SimModel = Simulink model name
% 2. Flight condition definition
% TrimParam.VelocitiesIni
% TrimParam.RatesIni
% TrimParam.AttitudeIni
% TrimParam.PositionIni
% TrimParam.FuelIni
% TrimParam.EngineSpeedIni
% TrimParam.Airspeed
% TrimParam.Altitude
% TrimParam.BankAngle
% TrimParam.Elevator
% TrimParam.Aileron
% TrimParam.Rudder
% TrimParam.Throttle
% TrimParam.Flap
% TrimParam.Mix
% TrimParam.Ign
% TrimParam.Winds
% 3. Miscellaneous parameters:
% TrimParam.StateIdx = state order index (order of states in Simulink block diagram usually
%
different than the desired state order - velocities, angular rates,
%
attitude angles, position, fuel, and engine speed)
% TrimParam.Options = trim function options
% TrimParam.SimOptions = the sim function options
% TrimParam.NAircraftStates = the size of the aircraft state vector,
%
currently 14
145
% TrimParam.NSimulinkStates = the size of the state vector in the Simulink
%
model (could be larger than aircraft state vector)
%
% The trim output structure has the following components:
% TrimOutput.States = Simulink model states at trim condition
% TrimOutput.Inputs = inputs at trim condition
% TrimOutput.Outputs = outputs at trim condition
% TrimOutput.Derivatives = model derivatives at trim condition
clear all;
clc;
fprintf('\nSetting initial trim parameters...');
% Simulation time settings
TrimParam.SampleTime = 0.04;
TrimParam.FinalTime = 60;
% Actuators
TrimParam.Flap = 0;
TrimParam.Mix = 13;
TrimParam.Ign = 1;
% Wind velocities
TrimParam.Winds = [0 0 0];
% Simulink model to trim
TrimParam.SimModel = 'aerosonde_trim';
fprintf('\nThe Simulink model %s.mdl will be trimmed.', TrimParam.SimModel);
% Get the sim options structure
TrimParam.SimOptions = simget(TrimParam.SimModel);
% Set the model inputs
TrimInput = [0 0 0 0.4];
%%% IDENTIFY THE ORDER OF AIRCRAFT STATES %%%
fprintf('\nIdentifying the order of the Simulink model states...');
146
% Set some easily-identifiable initial conditions
TrimParam.VelocitiesIni = [24 2 -1]';
TrimParam.RatesIni = [0.03 0.02 0.01]';
TrimParam.AttitudeIni = [0.1 0.05 0.6]';
TrimParam.PositionIni = [45*pi/180 -122*pi/180 1342]';
TrimParam.FuelIni = 2.567;
TrimParam.EngineSpeedIni = 5000*pi/30;
% Run Simulink model for a single sample period
[SimTime, SimStates, SimOutputs] = sim(TrimParam.SimModel, [0 TrimParam.SampleTime],
TrimParam.SimOptions, ...
[0 TrimInput; TrimParam.SampleTime TrimInput]);
% Find the state order
StateIni = [TrimParam.VelocitiesIni; TrimParam.RatesIni; TrimParam.AttitudeIni;
TrimParam.PositionIni; TrimParam.FuelIni; TrimParam.EngineSpeedIni];
TrimParam.NAircraftStates = length(StateIni);
TrimParam.NSimulinkStates = length(SimStates(1,:));
TrimParam.StateIdx = zeros(TrimParam.NAircraftStates,1);
for i=1:TrimParam.NAircraftStates
j = 1;
while TrimParam.StateIdx(i) == 0
value = SimStates(1,j);
if value == StateIni(i)
TrimParam.StateIdx(i) = j;
end
j = j + 1;
end
end
clear SimTime SimStates SimOutputs;
fprintf('done.');
%%% DETERMINE THE INITIAL GUESS FOR AIRCRAFT CONTROLS %%%
% Flight condition
fprintf('\n');
fprintf('\n Choose flight condition:');
fprintf('\n--------------------------\n');
TrimParam.Airspeed = input('Trim airspeed [m/s]: ');
TrimParam.Altitude = input('Trim altitude [m]: ');
147
TrimParam.BankAngle = input('Trim bank angle [rad]: ');
TrimParam.FuelIni = input('Fuel mass [kg]: ');
TrimParam.Flap = input('Flap setting [frac]: ');
% The initial conditions for the flight condition above
TrimParam.VelocitiesIni = [TrimParam.Airspeed 0 0]';
TrimParam.RatesIni = [0 0 0]';
TrimParam.AttitudeIni = [TrimParam.BankAngle 0 0]';
TrimParam.PositionIni = [45*pi/180 -122*pi/180 TrimParam.Altitude]';
TrimParam.EngineSpeedIni = 5000*pi/30;
% The trim error threshold
MaxErrTAS = 0.5;
MaxErrAlt = 4;
MaxErrBank = 0.1*pi/180;
% The control surface gains
KElevator = -0.02;
KAileron = 0.02;
KThrottle = -0.02;
fprintf('\nComputing the initial estimates for the trim inputs...');
GoodGuess = 0; Niter = 1;
while (~GoodGuess)&(Niter<30)
% Run Simulink model for a short time (10 s)
[SimTime, SimStates, SimOutputs] = sim(TrimParam.SimModel, [0 10], TrimParam.SimOptions, ...
[0 TrimInput; 10 TrimInput]);
% Compute errors in trim
ErrTAS = SimOutputs(end,1) - TrimParam.Airspeed;
ErrAlt = SimOutputs(end,7) - TrimParam.Altitude;
ErrBank = SimOutputs(end,4) - TrimParam.BankAngle;
fprintf('\nIteration #%2d, Airsp err = %6.2f m/s, Alt err = %8.2f m, phi err = %6.2f deg.', Niter,
ErrTAS, ErrAlt, ErrBank*180/pi);
% If all errors are within threshold
if (abs(ErrTAS)<MaxErrTAS)&(abs(ErrAlt)<MaxErrAlt)&(abs(ErrBank)<MaxErrBank)
% We are done with the initial guess
GoodGuess = 1;
else
148
% Adjust aircraft controls
TrimInput(1) = TrimInput(1) + KElevator * ErrTAS;
TrimInput(2) = TrimInput(2) + KAileron * ErrBank;
TrimInput(4) = TrimInput(4) + KThrottle * ErrAlt;
end
Niter = Niter + 1;
end
% Save initial guess
TrimParam.VelocitiesIni = SimStates(end,TrimParam.StateIdx(1:3))';
TrimParam.AttitudeIni = SimStates(end, TrimParam.StateIdx(7:9))';
TrimParam.EngineSpeedIni = SimStates(end, TrimParam.StateIdx(14));
TrimParam.Elevator = TrimInput(1);
TrimParam.Aileron = TrimInput(2);
TrimParam.Rudder = TrimInput(3);
TrimParam.Throttle = TrimInput(4);
clear SimTime SimStates SimOutputs TrimInput;
fprintf('\nDone. Initial guesses for trim inputs are:');
fprintf('\n Elevator = %6.4f', TrimParam.Elevator);
fprintf('\n Aileron = %6.4f', TrimParam.Aileron);
fprintf('\n Rudder = %6.4f', TrimParam.Rudder);
fprintf('\n Throttle = %6.4f', TrimParam.Throttle);
%%% PERFORM AIRCRAFT TRIM %%%
fprintf('\n');
fprintf('\nPerforming the aircraft trim...\n');
% Set initial guesses
StateIni = zeros(TrimParam.NSimulinkStates,1);
StateIni(TrimParam.StateIdx(1)) = TrimParam.VelocitiesIni(1);
StateIni(TrimParam.StateIdx(2)) = TrimParam.VelocitiesIni(2);
StateIni(TrimParam.StateIdx(3)) = TrimParam.VelocitiesIni(3);
StateIni(TrimParam.StateIdx(4)) = TrimParam.RatesIni(1);
StateIni(TrimParam.StateIdx(5)) = TrimParam.RatesIni(2);
StateIni(TrimParam.StateIdx(6)) = TrimParam.RatesIni(3);
StateIni(TrimParam.StateIdx(7)) = TrimParam.BankAngle;
StateIni(TrimParam.StateIdx(8)) = TrimParam.AttitudeIni(2);
StateIni(TrimParam.StateIdx(9)) = TrimParam.AttitudeIni(3);
StateIni(TrimParam.StateIdx(10)) = TrimParam.PositionIni(1);
StateIni(TrimParam.StateIdx(11)) = TrimParam.PositionIni(2);
StateIni(TrimParam.StateIdx(12)) = TrimParam.Altitude;
149
StateIni(TrimParam.StateIdx(13)) = TrimParam.FuelIni;
StateIni(TrimParam.StateIdx(14)) = TrimParam.EngineSpeedIni;
InputIni = [TrimParam.Elevator; TrimParam.Aileron; TrimParam.Rudder; TrimParam.Throttle];
OutputIni = [TrimParam.Airspeed; 0; 0; TrimParam.BankAngle; 0; 0; TrimParam.Altitude];
DerivIni = zeros(TrimParam.NSimulinkStates,1);
% Set indices of fixed parameters
StateFixIdx = [TrimParam.StateIdx(4) TrimParam.StateIdx(5) TrimParam.StateIdx(6)
TrimParam.StateIdx(7) TrimParam.StateIdx(12) TrimParam.StateIdx(13)];
InputFixIdx = [];
OutputFixIdx = [1 2 4 7];
DerivFixIdx = [TrimParam.StateIdx(1) TrimParam.StateIdx(2) TrimParam.StateIdx(3)
TrimParam.StateIdx(4) TrimParam.StateIdx(5) TrimParam.StateIdx(6) ...
TrimParam.StateIdx(7) TrimParam.StateIdx(8) TrimParam.StateIdx(9) TrimParam.StateIdx(12)
TrimParam.StateIdx(14)];
% Set optimization parameters
TrimParam.Options(1) = 1;
% show some output
TrimParam.Options(2) = 1e-6; % tolerance in X
TrimParam.Options(3) = 1e-6; % tolerance in F
TrimParam.Options(4) = 1e-6;
TrimParam.Options(14) = 5000; % max iterations
% Trim the airplane
[TrimOutput.States,TrimOutput.Inputs,TrimOutput.Outputs,TrimOutput.Derivatives] =
trim(TrimParam.SimModel,StateIni,InputIni,OutputIni,...
StateFixIdx,InputFixIdx,OutputFixIdx,DerivIni,DerivFixIdx,TrimParam.Options);
% Print the trim results
fprintf('\nFinished. The trim results are:');
fprintf('\nINPUTS:');
fprintf('\n Elevator = %6.4f', TrimOutput.Inputs(1));
fprintf('\n Aileron = %6.4f', TrimOutput.Inputs(2));
fprintf('\n Rudder = %6.4f', TrimOutput.Inputs(3));
fprintf('\n Throttle = %6.4f', TrimOutput.Inputs(4));
fprintf('\nSTATES:');
fprintf('\n u = %6.2f m/s', TrimOutput.States(TrimParam.StateIdx(1)));
fprintf('\n v = %6.2f m/s', TrimOutput.States(TrimParam.StateIdx(2)));
fprintf('\n w = %6.2f m/s', TrimOutput.States(TrimParam.StateIdx(3)));
150
fprintf('\n p = %6.2f deg/s', TrimOutput.States(TrimParam.StateIdx(4))*180/pi);
fprintf('\n q = %6.2f deg/s', TrimOutput.States(TrimParam.StateIdx(5))*180/pi);
fprintf('\n r = %6.2f deg/s', TrimOutput.States(TrimParam.StateIdx(6))*180/pi);
fprintf('\n phi = %6.2f deg', TrimOutput.States(TrimParam.StateIdx(7))*180/pi);
fprintf('\n theta = %6.2f deg', TrimOutput.States(TrimParam.StateIdx(8))*180/pi);
fprintf('\n psi = %6.2f deg', TrimOutput.States(TrimParam.StateIdx(9))*180/pi);
% Geographic position is irrelevant
%fprintf('\n Lat = %8.4f deg', TrimOutput.States(TrimParam.StateIdx(10))*180/pi);
%fprintf('\n Lon = %8.4f deg', TrimOutput.States(TrimParam.StateIdx(11))*180/pi);
fprintf('\n Alt = %6.2f m', TrimOutput.States(TrimParam.StateIdx(12)));
fprintf('\n Fuel = %6.2f kg', TrimOutput.States(TrimParam.StateIdx(13)));
fprintf('\n Engine = %6.0f rot/min', TrimOutput.States(TrimParam.StateIdx(14))*30/pi);
fprintf('\nOUTPUTS:');
fprintf('\n Airspeed = %6.2f m/s', TrimOutput.Outputs(1));
fprintf('\n Sideslip = %6.2f deg', TrimOutput.Outputs(2)*180/pi);
fprintf('\n AOA = %6.2f deg', TrimOutput.Outputs(3)*180/pi);
fprintf('\n Bank = %6.2f deg', TrimOutput.Outputs(4)*180/pi);
fprintf('\n Pitch = %6.2f deg', TrimOutput.Outputs(5)*180/pi);
fprintf('\n Heading = %6.2f deg', TrimOutput.Outputs(6)*180/pi);
fprintf('\n Altitude = %8.2f m', TrimOutput.Outputs(7));
% Extract the linear model
fprintf('\n \nExtracting aircraft linear model...\n');
% Perturbation level
LinParam(1) = 10^-8;
[A, B, C, D] = linmod(TrimParam.SimModel, TrimOutput.States, TrimOutput.Inputs, LinParam);
% Longitudinal dynamics
% States: u w q theta h omega
% Inputs: elevator throttle
% Outputs: Va alpha q theta h
Alon = [
A(TrimParam.StateIdx(1), [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5)
TrimParam.StateIdx(8) TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
A(TrimParam.StateIdx(3), [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5)
TrimParam.StateIdx(8) TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
A(TrimParam.StateIdx(5), [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5)
TrimParam.StateIdx(8) TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
A(TrimParam.StateIdx(8), [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5)
TrimParam.StateIdx(8) TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
151
A(TrimParam.StateIdx(12), [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5)
TrimParam.StateIdx(8) TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
A(TrimParam.StateIdx(14), [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5)
TrimParam.StateIdx(8) TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
];
Blon = [
B(TrimParam.StateIdx(1), [1 4])
B(TrimParam.StateIdx(3), [1 4])
B(TrimParam.StateIdx(5), [1 4])
B(TrimParam.StateIdx(8), [1 4])
B(TrimParam.StateIdx(12), [1 4])
B(TrimParam.StateIdx(14), [1 4])
];
Clon = [
C(1, [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5) TrimParam.StateIdx(8)
TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
C(3, [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5) TrimParam.StateIdx(8)
TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
zeros(2, 6)
C(7, [TrimParam.StateIdx(1) TrimParam.StateIdx(3) TrimParam.StateIdx(5) TrimParam.StateIdx(8)
TrimParam.StateIdx(12) TrimParam.StateIdx(14)])
];
Clon(3,3) = 1; Clon(4,4) = 1;
fprintf('\n');
fprintf('\n Longitudinal Dynamics');
fprintf('\n-----------------------');
fprintf('\n State vector: x = [u w q theta h Omega]');
fprintf('\n Input vector: u = [elevator throttle]');
fprintf('\n Output vector: y = [Va alpha q theta h]');
fprintf('\n State matrix: A = \n');
disp(Alon);
fprintf('\n Control matrix: B = \n');
disp(Blon);
fprintf('\n Observation matrix: C = \n');
disp(Clon);
% Eigenvalue analysis
eiglon = eig(Alon);
for i=1:length(eiglon)
152
if imag(eiglon(i))>0
[wd, T, wn, zeta] = eigparam(eiglon(i));
fprintf('\n Eigenvalue: %6.4f +/- %6.4f i', real(eiglon(i)), imag(eiglon(i)));
fprintf('\n Damping = %6.4f, natural frequency = %6.4f rad/s, period = %8.4f s', zeta, wn, T);
elseif imag(eiglon(i))==0
fprintf('\n Eigenvalue: %6.4f', eiglon(i));
fprintf('\n Time constant = %6.4f s', -1/eiglon(i));
end
end
% Lateral-directional dynamics
% States: v p r phi psi
% Inputs: aileron rudder
% Outputs: beta p r phi psi
Alat = [
A(TrimParam.StateIdx(2), [TrimParam.StateIdx(2) TrimParam.StateIdx(4) TrimParam.StateIdx(6)
TrimParam.StateIdx(7) TrimParam.StateIdx(9)])
A(TrimParam.StateIdx(4), [TrimParam.StateIdx(2) TrimParam.StateIdx(4) TrimParam.StateIdx(6)
TrimParam.StateIdx(7) TrimParam.StateIdx(9)])
A(TrimParam.StateIdx(6), [TrimParam.StateIdx(2) TrimParam.StateIdx(4) TrimParam.StateIdx(6)
TrimParam.StateIdx(7) TrimParam.StateIdx(9)])
A(TrimParam.StateIdx(7), [TrimParam.StateIdx(2) TrimParam.StateIdx(4) TrimParam.StateIdx(6)
TrimParam.StateIdx(7) TrimParam.StateIdx(9)])
A(TrimParam.StateIdx(9), [TrimParam.StateIdx(2) TrimParam.StateIdx(4) TrimParam.StateIdx(6)
TrimParam.StateIdx(7) TrimParam.StateIdx(9)])
];
Blat = [
B(TrimParam.StateIdx(2), [2 3])
B(TrimParam.StateIdx(4), [2 3])
B(TrimParam.StateIdx(6), [2 3])
B(TrimParam.StateIdx(7), [2 3])
B(TrimParam.StateIdx(9), [2 3])
];
Clat = [
C(2, [TrimParam.StateIdx(2) TrimParam.StateIdx(4) TrimParam.StateIdx(6) TrimParam.StateIdx(7)
TrimParam.StateIdx(9)])
zeros(4, 5)
];
Clat(2,2) = 1; Clat(3,3) = 1;
Clat(4,4) = 1; Clat(5,5) = 1;
153
fprintf('\n');
fprintf('\n Lateral-directional Dynamics');
fprintf('\n------------------------------');
fprintf('\n State vector: x = [v p r phi psi]');
fprintf('\n Input vector: u = [aileron rudder]');
fprintf('\n Output vector: y = [beta p r phi psi]');
fprintf('\n State matrix: A = \n');
disp(Alat);
fprintf('\n Control matrix: B = \n');
disp(Blat);
fprintf('\n Observation matrix: C = \n');
disp(Clat);
% Eigenvalue analysis
eiglat = eig(Alat);
for i=1:length(eiglat)
if imag(eiglat(i))>0
[wd, T, wn, zeta] = eigparam(eiglat(i));
fprintf('\n Eigenvalue: %6.4f +/- %6.4f i', real(eiglat(i)), imag(eiglat(i)));
fprintf('\n Damping = %6.4f, natural frequency = %6.4f rad/s, period = %8.4f s', zeta, wn, T);
elseif (imag(eiglat(i))==0)&(real(eiglat(i))~=0)
fprintf('\n Eigenvalue: %6.4f', eiglat(i));
fprintf('\n Time constant = %6.4f s', -1/eiglat(i));
end
end
154