Trimagem do sistema não linear do F16

Após a transcrição dos dados presentes no livro do Lewis, 3ª edição, relativos ao F16, os coeficientes aerodinâmicos simplificados estão presentes nas Figuras 1 a 4. Nelas é possível verificar a existência de superfícies suaves (por exemplo o coeficiente da Figura 1.b) e em outros casos a superfície em si apresenta pontos acentuados de descontinuidade (majoritariamente nos coeficientes expostos nas Figuras 3 e 4).

Para confrontar os resultados das simulações linear e não linear, é importante então que inicialmente seja definido o ponto de operação desejado. Esse ponto de operação desejado será tomado como referência para o algoritmo de trimagem. Nesse ponto, será encontrado o equilíbrio teórico, mesmo que momentâneo/instantâneo, nas situações em que se supõe a não alteração de massa por consumo de combustível e/ou momentos de inércia (seja por variação de massa e/ou deflexão de empenas de atuação do sistema).

Figura 1: (a) Coeficientes de Momento de Rolamento e (b) de Pitching.

Figura 2: (a) Coeficientes de Momento de Yawing e de Força no Eixo X.

Figura 3: (a) Coeficientes de Momento de Rolamento devido ao Rudder (b) e do Aileron.

Figura 4: (a) Coeficientes de Momento de Yawing devido ao Aileron e (b) devido ao Rudder.

O algoritmo de trimagem utiliza a função fminsearch do Matlab para tal. É importante destacar que essa função em específico realiza a minimização do sistema sem considerar qualquer tipo de restrições nas variáveis. A depender de como se comporta a semente (seed) de iniciação de busca, é provável que o algoritmo pare em alguma condição de mínimo local. Dessa forma é importante ao projetista que se conheça um pouco da dinâmica do sistema ao qual você deseja realizar a trimagem de forma a possibilitar uma semente inicial que seja condizente com uma factibilidade do seu sistema.

Em atualização recente, eu utilizei a função fmincon, onde que insiro algumas restrições referentes à parte de range de atuação dos ailerons, elevator e rudder. E nesse caso, o sistema tem sido forçado a encontrar um mínimo global e não apenas um local (não precisaria da ação de um agente externo ou de um grande conhecedor da dinâmica da aeronave para vislumbrar a inicialização das variáveis do sistema).

O programa para trimagem está descrito a partir da Figura 5 e explicado passo a passo o seu funcionamento. Inicialmente todos os arquivos que dizem respeito ao F-16 são carregados (vide a carga das pastas e a adição das mesmas no caminho – path – salvo). Após a carga das pastas e arquivos, surge para o programador a tela de decisão (vide Figura 6) para a condição de voo (no caso está implementado apenas a condição de "reto nivelado" e "curva coordenada"). Tal tela é possibilitada pelas linhas da Figura 7. Após essa seleção, o programa em si parte para um switch-case a depender da situação desejada pelo programador.

Figura 5: Carrega as pastas e arquivos que serão utilizados no programa de simulação.

Figura 6: Seleção da condição de voo. Foi comentada a situação de rolamento estacionário para desenvolvimento futuro.

Uma vez na função switch-case, o programa segue a rotina da Figura 7.

Figura 7: Caracterização das variáveis a serem otimizadas pela função fmincon, limites superiores e inferiores. Observe o case para o throtle, que foi colocado um parâmetro ed 1.5 de limite superior (isso considera o afterburn).

Além disso, se o seu avião está voando em reto nivelado, é interessante também que haja (ou seja, não exista ângulo entre o vetor velocidade total e o plano NED). Além disso, não há nenhum nível de rolamento do seu avião, assim . O que sobra em termos de variáveis que possibilitam a otimização do sistema dizem respeito à , h e (as outras variáveis são determinadas pela restrição direta do problema). Após essa iniciação pelo usuário, é necessário que o algoritmo inicie a otimização das variáveis. No caso, dado as restrições físicas do próprio problema, as variáveis livres são dadas por:

O algoritmo, então, a partir dessas variáveis (chute inicial e em seguida realimentação dos valores previstos) põe tais valores nas equações de estado e calcula função de custo. O algoritmo checa se houve ou não convergência dos valores para algum ponto de mínimo e fica nessa busca até encontrar algum de mínimo. A variável

será utiliza dentro da função fnCusto (a qual será detalhada a seguir) para travar algumas restrições físicas do sistema. Além disso, é importante que se dê um chute inicial nas variáveis livres.

Figura 8: Função fnCusto.

A função fnCusto recebe como parâmetros o vetor s e o varargin. Eu inicialmente quebro esse vetor s nas variáveis às quais eu desejo minimizar, pois tais variáveis entrarão na função que descreve a dinâmica do sistema do avião (e que aceita em específico o vetor de estados x, o vetor de controle u e o tempo t). Equivalentemente, faço a mesma coisa para o vetor varargin. Em termos de restrições existentes para o estado, o algoritmo identifica qual a situação de voo desejada (dado pela variável Indice) e já deixa explícitas as restrições do sistema. Após deixado o problema bem amarrado em termos de restrições, o algoritmo calcula a derivada do vetor de estados dadas as condições do voo escolhida e as restrições físicas existentes. Em seguida, calcula-se o custo quadrático do sistema em cima dos resultados obtidos para a derivada do vetor de estados .

Após a realização da otimização, coleto os valores obtidos .

Figura 9: Pego nas variáveis de estado e de controle os valores que minimizam a função fnCusto.

Após a obtenção das variáveis minimizantes, é necessário que se calcule também as outras variáveis de estado (no caso aquelas que estão obviamente restritas). Esse passo de recálculo, dados as restrições, é realizado na Figura 10. Nela, a depender da condição de voo, realiza-se a mesma toada de equações de restrições de variáveis existentes dentro da função fnCusto. Após a realização do cálculo do sistema, eu calculo a variável de estado x(13) e salvo os resultados em um arquivo trimados.mat.

Exemplo de aplicação (o mesmo que está no livro...): com a finalidade de elucidar o procedimento de trimagem para uma determinada condição de voo, eu irei considerar a seguinte situação de voo

• 
• 
• 
• 
• Voo reto nilevado, o que implica em .

Ou seja, nada mais é do que a condição de voo 1 do nosso problema.

Figura 10: Calculo as variáveis de estado restantes que minimizam a função fnCusto de acordo com as restrições de voo.

Para a situação de trimagem desejada, o resultado obtido na otimização encontra-se apresentado na próxima figura. É importante deixar claro que as 9 primeiras variáveis são as variáveis de estados (a parte de posição para navegação não é importante para este caso de trimagem) e os 4 últimos são referentes aos parâmetros de controle.

Figura 11: Vetor de estados e controles trimados. É possível observar o custo resultante da otimização.

Esse procedimento pode ser ser repetido para diversas condições de voo para dentro do envelope de voo (altitude, velocidade, condição de variação de massa) a ser considerado para o seu projeto de aeronave e/ou drone.

Alguns pontos que deve ser levados em consideração, os valores de trimagem para "AIL" e RDR" são muito pequenos e podem ser considerados como nulos na trimagem. Isso fatalmente irá ocorrer em outras situações de trimagem do seu sistema. Fique a vontade, por exemplo, de zerá-los quando da determianção do seus parâmetros linearizados. Para os desenvolvimentos que irei fazer, eu vou tomar esses valores obtidos diretamente na trimagem, mesmo que muito pequeno, quando na obtenção das matrizes de linearização do meu sistema para desenvolver malhas de controle.

É interessante observar que esse procedimento pode ser adotado de maneira similar quando se tem disponível um modelo não linear e se deseja encontrar algum ponto de trimagem para projeto ou ajustes em malhas de controle.

Modelo não linear F16

O modelo não-linear encontra-se disponível no livro do Stevens e Lewis, Airfract Control and Simulation. Possuo as 2ª e 3ª edições do livro, e não há variação no código existente entre elas. O código disponibilizado no livro para a dinâmica não linear do F16 está escrito em Fortran e para facilitar um pouco a vida eu a transcrevi para o Matlab. Hoje em dia com a IA você consegue fazer isso num tapa, mas a versão por mim programada pode ser encontrada no meu github por meio deste link. O modelo utiliza o referencial NED (North-East-Down). Você pode, por exemplo, utilizar o ENU (East-North-Up), caso deseje, por questões de opção de referencial, mas isso não consta no programação realizada em Matlab.

Há de convir que eu procurei tentar separar o que pega no sistema de modelagem matemática. Tentei separar em subsistemas de modo a deixar a programação separada em massa, momentos de inercia, coeficientes aerodinâmicos, influência de momento rotor (turbina, por exemplo), entre outros, e deixar as equações de força, momento, cinemática e dinâmica um pouco mais direto para quem deseja entender o passo a passo da modelagem.

A modelagem considera a terra plana, mas dá para realizar algumas adaptações nesse código de modo a considerar a terra esférica (e as influências de Corioli e força centrífuga) na modelagem. Isso é um preciosismo, mas para a tarefa de obtenção de malhas de controle por agora não tem muuuuito ganho.

Considerações de carater geral:

• As condicionantes de massa em uma aeronave, drone com propulsão sólida/líquida varia com o passar do tempo, pelo consumo do agente propulsor. No caso da dinâmica F16 modelado e programado no livro do Lewis, essa massa não varia. No seu programa mais refinado, para a dinâmica que você irá desenvolver de modo mais profissional, leve em consideração as variações de massa total do seu sistema. Se seu sistema possui baterias elétricas, essa variação não tem fundamentação, pois não há variação mássica e apenas consumo de energia;
• Equivalentemente, há alteração de momento de inércia do seu sistema com variação da condição de voo (flaps baixados, variações devido às deflexões de empenas, devido à própria variação de massa, alijamento de algum partes não necessárias mais durante o voo). Isso deve ser levado em consideração na sua abordagem quando da simulação para o caso mais geral;
• Além disso, em se considerando a terra esférica, deve-se considerar variações no parâmetro de gravidade com a latitude, como mesmo abordado no livro do Lewis;
• Os parâmetros aerodinâmicos (tabelas aerodinâmicas) e as de desempenho do turbina do F16 vêm diretamente das tabelas presentes no livro do Lewis. No caso de sua simulação de aeronave, drone você tem que derivar essas tabelas aerodinâmicas por CFD, análise numérica ou qualquer outro tipo aplicável. Aqui eu não fiz juízo de valor sobre se elas estão corretas ou se os seus parâmetros possuem equivalência com o mundo real. Apenas as usei;
• A equação para o F16 já foi construída tendo por base o parâmetro t, que num caso geral toma por referência o consumo ou alguma máquina de estados do sistema em condição após certo instante de tempo, por exemplo.

Existem alguns códigos disponíveis na internet, até com novas tabelas aerodinâmicas. Existe também alguns trabalhos, livros e afins que consideram o modelo do F16 uma base para desenvolvimento e estudos na área de controle, otimização numérica, simulação de dinâmica etc. Esse é o código, como dito anteriormente, escrito por mim, e se tiver erro, faz parte. Me aponte onde, por favor, que eu ajusto tão logo ele seja a mim chegado.

O modelo foi escrito basicamente seguindo as diretrizes disponíveis no livro. O código a seguir irá esclarecer as etapas do desenvolvimento. É possível observar que existem alguns comentários no código, que fazem referência direta à equação tal do livro, disponível na página tal. Segue o código.

Figura 1: Num contexto mais geral, a sua função matemática irá cuspir uma massa, peso e matriz J, referente ao momento de inércia do sistema, condizente com a condição de voo em que o sistema se encontra. No fim das contas, é um parâmetro dependente do parâmetro t.

Figura 2: Parâmetros específicos da aeronave F16 apresentados no livro do Lewis. Num sistema complexo, esses parâmetros podem, e com certeza, devem variar.

Figura 3: Estado instantâneo do sistema para as variáveis de estado e controle.

Figura 4: Análise do sistema num pré-cálculo de variáveis que serão empregadas no transcorrer do código.

Figura 5: Modelo aerodinâmico do sistema na condição ISA + 0 apresentado no livro do Lewis. Para um modelo mais completo, é interessante realizar a inserção do parâmetro de temperatura local para adaptar para as demais ISA's.

Figura 6: Modelo propulsivo do F16 apresentado pelo Lewis. Aqui o buraco é mais embaixo num sistema complexo, em se considerando rendimento, momento angular induzido entre outros parâmetros. No livro do Lewis está tudo basicamente em formato de tabela.

Figura 7: Etapa de processamento dos parâmetros aerodinâmicos. Em sistemas não lineares é comum encontrar tabelas e análise nesse formato, mais ou menos complexas, e levando em considerações as diversas fases de uma aeronave, drone em seus mais diversos estágios de voo.

Figura 8: Equações de Força do sistema.Ela está diferente daquela apresentada no livro do Lewis de modo a deixar o desenvolvimento mais direto com as equações disponíveis.

Figura 9: Equações de cinemática. Aqui a coisa não muda muito do que já existe. Tem que se levar em considerações a possibilidade de divisão por zero em alguns termos. No futuro, interessante modificar para quaternions para evitar problemas numéricos.

Figura 10: Equações de momento. Observe a consideração de momento induzida pelos motores. No caso mais geral, há também a dinâmica do momento de inércia , que será derivado a posteriori em outro post com mais detalhes de como o fazer matematicamente.

Figura 11: Equações de navegação, tendo por referência o padrão NED. Fica a critério do cliente aqui apresentar os parâmetros derivados.

Figura 12: Saídas de interesse, a critério do desenvolvedor e das malhas as quais ele deseja obter para o desempenho do seu sistema.

Bom, esse é omodelo simplificado do F16 disponível no livro do Lewis. É possível, por exemplo, criar um s-function no Simulink ou mesmo derivar, bloco a bloco, para no próprio Simulink, para deixar a coisa mais impessoal. O código gera os parâmetros apresentados no livro, quando da verificação de modelo. Há um pequeno "erro" para a variável , desprezíveis quando da análise de um contexto mais geral e de testes para aprendizado do problema. Vi alguns códigos e coisas do tipo em C++ na internet e também em Matlab. Quero deixar disponível no blog apenas para consulta futura minha, caso se faça necessário.

Ah, para entender melhor as funções que são chamadas durante o programa, clica no link, deixa de preguiça e olha dentro das pastas para você entender o que foi feito. É isso aí. Vlw.