INSTITUTO TECNOLÓGICO DE AERONÁUTICA

Tarik Carvalho Bussab

Simulação de um Sistema de Controle

de um Quadricóptero

Trabalho de Graduação

2014

Mecânica-Aeronáutica

CDU 629.73.017.2

TARIK CARVALHO BUSSAB

SIMULAÇÃO DE UM SISTEMA DE CONTROLE

DE UM QUADRICÓPTERO

Orientador

Prof. Dr. Davi Antônio dos Santos. (ITA)

Divisão de Engenharia Mecânica-Aeronáutica

SÃO JOSÉ DOS CAMPOS

INSTITUTO TECNOLÓGICO DE AERONÁUTICA

2014

Dados Internacionais de Catalogação-na-Publicação (CIP)

Divisão de Informação e Documentação

Bussab, Tarik Carvalho

Simulação de um Sistema de Controle de um Quadricóptero / Tarik Carvalho Bussab

São José dos Campos, 2014.

57f.

Trabalho de Graduação – Divisão de Engenharia Mecânica-Aeronáutica –

Instituto Tecnológico de Aeronáutica, 2014. Orientadores: Prof. Dr. Davi Antônio dos Santos.

1. Quadricóptero. 2. Controle. 3. Simulação. II. Instituto Tecnológico de Aeronáutica. III.

Simulação de um Sistema de Controle de um Quadricóptero.

REFERÊNCIA BIBLIOGRÁFICA

BUSSAB, Tarik Carvalho. Simulação de um Sistema de Controle de um Quadricóptero. 2014.

57 folhas. Trabalho de Conclusão de Curso. (Graduação) – Instituto Tecnológico de

Aeronáutica, São José dos Campos.

CESSÃO DE DIREITOS

NOME DO AUTOR: Tarik Carvalho Bussab

TÍTULO DO TRABALHO: Simulação de um Sistema de Controle de um Quadricóptero

TIPO DO TRABALHO/ANO: Graduação / 2014

É concedida ao Instituto Tecnológico de Aeronáutica permissão para reproduzir cópias

deste trabalho de graduação e para emprestar ou vender cópias somente para propósitos

acadêmicos e científicos. O autor reserva outros direitos de publicação e nenhuma parte

desta monografia de graduação pode ser reproduzida sem a autorização do autor.

Tarik Carvalho Bussab

Rua Visconde de Taunay, 443, ap. 54

São Paulo – SP, 04726-010

Dedico este trabalho à minha família e

aos meus professores, que me deram

base para realizar este trabalho.

Agradecimentos

Agradeço primeiramente aos meus pais, que me ensinaram os valores do estudo, da

dedicação e da força.

Foi fundamental a orientação do Prof. Dr. Davi Antônio dos Santos a este trabalho. Sua

maneira de ensinar aos poucos e de ensinar o aluno a pensar é muito valiosa.

Agradeço também aos Prof. Ijar Milagre da Fonseca e Prof. Rafael Thiago Luiz Ferreira,

que deram contribuições e comentários valiosos para a realização deste trabalho.

“The gods envy us. They envy us because we are mortal,

because any moment may be our last. Everything is more beautiful because we are

doomed. You will never be lovelier than you are now.

We will never be here again.”

Resumo

Este trabalho se dedicará à modelagem, simulação e controle de atitude e posição de um robô aéreo

do tipo quadricóptero. É bem conhecido que o veículo em questão apresenta uma dinâmica

subatuada, pois contém seis graus de libertade (três de rotação e três de translação), enquanto está

sujeito a apenas quatro variáveis de atuação independentes (três componentes de torque e um

componente de força perpendicular ao plano dos rotores). Devido a essa característica, o

quadricóptero não poderia ser controlado independentemente em seus seis graus de liberdade. Sem

embargo, isso não consiste numa dificuldade, desde que o interesse central é normalmente

controlar independentemente apenas quatro graus de liberdade: os três componentes de posição e o

ângulo de heading. O sistema que realiza tal controle é composto de duas malhas aninhadas, em

que a interna se ocupa do controle de atitude e a externa é responsável por controlar a posição. A

modelagem dinâmica será realizada mediante a formulação Newton-Euler. Na malha interna (de

atitude) far-se-á uso de três controladores proporcionais-derivativos (PD) desacoplados. O controle

da malha externa (de posição) será realizado pela lei de controle proposta em [3], que consiste num

controlador PD saturado em combinação com uma linearização por realimentação de estados. O

sistema descrito será simulado no ambiente MATLAB/Simulink. A simulação consiste em o robô

ser comandado a passar por dois pontos no espaço e depois retornar à sua posição inicial.

Palavras-chave: Quadricóptero, Simulação, Controlador de Atitude, Controlador de

Posição, Dinâmica, Modelagem

Summary

This work was dedicated to the modeling, simulation and attitude and position control of a aerial

vehicle of type quadcopter. It is well known that such vehicle has a under controlled dynamics,

since it has six degrees of freedom (three of rotation and three of translation), but it is only subject

to four independent control variables (three torque components and one component of the force

normal to the plane that contains the rotors). Due to this characteristic, the quadcopter cannot be

controlled independently in it six degrees of freedom. Nevertheless, this is not an issue since the

main interest is generally to independently control only four degrees of freedom: three position

components and the heading angle. The system that undertakes this control comprises two nested

loops, where the inner one carries attitude control and the outer one is responsible for the position

control. The dynamic modeling will be derived via Newton-Euler’s formulations. In the inner loop

(of attitude) three uncoupled proportional-derivative (PD) controllers will be implemented. The

outer loop control (of position) will be based on the control laws proposed in [3], which consists of

a saturated PD controller together with a linearization by the states feedback. The system described

will be simulated on the MATLAB/Simulink environment. The simulation consists in the robot

being commanded to visit two points in space and then return to its initial position.

Keywords: Quadcopter, Simulation, Attitude Controller, Position Controller, Dynamics,

Modeling.

Lista de ilustrações

Figura 3.1 – Esquema do quadricóptero e definição dos sistemas de coordenas..............................17

Figura 3.2 – Rotação do sistema de coordenadas segundo Eixo e Ângulo de Euler.........................20

Figura 4.1 – Desenho do quadricóptero em CAD.............................................................................24

Figura 5.1 – Ambiente Simulink para testes da dinâmica.................................................................26

Figura 5.2 – Pairar | 𝑥𝐶𝑀....................................................................................................................28

Figura 5.3 – Pairar | 𝑦𝐶𝑀....................................................................................................................28

Figura 5.4 – Pairar | 𝑧𝐶𝑀....................................................................................................................28

Figura 5.5 – Pairar | �̇�𝐶𝑀....................................................................................................................28

Figura 5.6 – Pairar | �̇�𝐶𝑀....................................................................................................................28

Figura 5.7 – Pairar | �̇�𝐶𝑀....................................................................................................................28

Figura 5.8 – Pairar | 𝜙........................................................................................................................28

Figura 5.9 – Pairar | 𝜃........................................................................................................................28

Figura 5.10 – Pairar | 𝜓......................................................................................................................28

Figura 5.11 – Pairar | 𝜔𝑥....................................................................................................................29

Figura 5.12 – Pairar | 𝜔𝑦....................................................................................................................29

Figura 5.13 – Pairar | 𝜔𝑧....................................................................................................................29

Figura 5.14 – Subir | 𝑥𝐶𝑀...................................................................................................................29

Figura 5.15 – Subir | 𝑦𝐶𝑀...................................................................................................................29

Figura 5.16 – Subir | 𝑧𝐶𝑀...................................................................................................................29

Figura 5.17 – Subir | �̇�𝐶𝑀...................................................................................................................30

Figura 5.18 – Subir | �̇�𝐶𝑀...................................................................................................................30

Figura 5.19 – Subir | �̇�𝐶𝑀...................................................................................................................30

Figura 5.20 – Subir | 𝜙.......................................................................................................................30

Figura 5.21 – Subir | 𝜃.......................................................................................................................30

Figura 5.22 – Subir | 𝜓......................................................................................................................30

Figura 5.23 – Subir | 𝜔𝑥.....................................................................................................................30

Figura 5.24 – Subir | 𝜔𝑦.....................................................................................................................30

Figura 5.25 – Subir | 𝜔𝑧.....................................................................................................................30

Figura 5.26 – Lateral | 𝑥𝐶𝑀................................................................................................................31

Figura 5.27 – Lateral | 𝑦𝐶𝑀................................................................................................................31

Figura 5.28 – Lateral | 𝑧𝐶𝑀................................................................................................................31

Figura 5.29 – Lateral | �̇�𝐶𝑀................................................................................................................32

Figura 5.30 – Lateral | �̇�𝐶𝑀................................................................................................................32

Figura 5.31 – Lateral | �̇�𝐶𝑀................................................................................................................32

Figura 5.32 – Lateral | 𝜙....................................................................................................................32

Figura 5.33 – Lateral | 𝜃....................................................................................................................32

Figura 5.34 – Lateral | 𝜓....................................................................................................................32

Figura 5.35 – Lateral | 𝜔𝑥..................................................................................................................32

Figura 5.36 – Lateral | 𝜔𝑦..................................................................................................................32

Figura 5.37 – Lateral | 𝜔𝑧..................................................................................................................32

Figura 5.38 – Yaw | 𝑥𝐶𝑀....................................................................................................................33

Figura 5.39 – Yaw | 𝑦𝐶𝑀....................................................................................................................33

Figura 5.40 – Yaw | 𝑧𝐶𝑀.....................................................................................................................33

Figura 5.41 – Yaw | �̇�𝐶𝑀....................................................................................................................33

Figura 5.42 – Yaw | �̇�𝐶𝑀....................................................................................................................33

Figura 5.43 – Yaw | �̇�𝐶𝑀.....................................................................................................................33

Figura 5.44 – Yaw | 𝜙........................................................................................................................34

Figura 5.45 – Yaw | 𝜃.........................................................................................................................34

Figura 5.46 – Yaw | 𝜓........................................................................................................................34

Figura 5.47 – Yaw | 𝜔𝑥......................................................................................................................34

Figura 5.48 – Yaw | 𝜔𝑦......................................................................................................................34

Figura 5.49 – Yaw | 𝜔𝑧.......................................................................................................................34

Figura 6.1 – Esquema simplificado das malhas fechadas do controlador do quadricóptero.............35

Figura 6.2 – Esquema auxiliar para determinação das leis de controle.............................................36

Figura 7.1 – Esquema simplificado das malhas fechadas do controlador do quadricóptero.............41

Figura 7.2 – Esquema completo das malhas fechadas do controlador do quadricóptero..................42

Figura 7.3 – Comando de 𝑥𝐶𝑀...........................................................................................................43

Figura 7.4 – Comando de 𝑦𝐶𝑀...........................................................................................................43

Figura 7.5 – Comando de 𝑧𝐶𝑀...........................................................................................................43

Figura 8.1 – Imagem em perspectiva da trajetória final....................................................................46

Figura 8.2 – Final | 𝑥𝐶𝑀.....................................................................................................................47

Figura 8.3 – Final | 𝑦𝐶𝑀.....................................................................................................................47

Figura 8.4 – Final | 𝑧𝐶𝑀.....................................................................................................................47

Figura 8.5 – Final | �̇�𝐶𝑀.....................................................................................................................47

Figura 8.6 – Final | �̇�𝐶𝑀.....................................................................................................................47

Figura 8.7 – Final | �̇�𝐶𝑀.....................................................................................................................47

Figura 8.8 – Final | 𝜙.........................................................................................................................47

Figura 8.9 – Final | 𝜃.........................................................................................................................47

Figura 8.10 – Final | 𝜓.......................................................................................................................47

Figura 8.11 – Final | 𝜔𝑥.....................................................................................................................48

Figura 8.12 – Final | 𝜔𝑦.....................................................................................................................48

Figura 8.13 – Final | 𝜔𝑧.....................................................................................................................48

Figura 8.14 – Ampliação da Fig. 8.2.................................................................................................48

Figura 8.15 – Ampliação da Fig. 8.3.................................................................................................48

Figura 8.16 – Ampliação da Fig. 8.4.................................................................................................48

Lista de abreviaturas e siglas

CM..........................................Centro de Massa

MRP........................................Parâmetros de Rodrigues Modificados

EDO........................................Equação Diferencial Ordinária

PD...........................................Proporcional-Derivativo

CAD.......................................Computer Aided Design

Lista de símbolos

𝑆𝐼..........................................Sistema de coordenadas fixo inercial

𝑆𝑅.........................................Sistema de coordenadas de referência

𝑆𝐶.........................................Sistema de coordenadas fixo no corpo

𝒓𝑪𝑴......................................Vetor posição do CM do quadricóptero escrito na base 𝑆𝐼

𝒇𝒊..........................................Força de propulsão gerada pela hélice 𝑖, com 𝑖 = 1, 2, 3, 4

𝑷..........................................Peso do corpo

ℓ...........................................Braço, a distância entre cada hélice e o CM

𝒎.........................................MRP representando a atitude de 𝑆𝐶 com respeito a 𝑆𝑅

𝝎..........................................Velocidade angular, escrita na base 𝑆𝐶, de 𝑆𝐶 com relação a 𝑆𝑅

�̂�...........................................Eixo unitário de rotação segundo a representação de atitude por Eixo e

Ângulo de Euler

𝜑...........................................Ângulo de rotação ao redor de �̂� segundo a representação de atitude

por Eixo e Ângulo de Euler

𝐼3..........................................Matriz Identidade de terceira ordem

[𝑚×].....................................Matriz que equivale ao produto vetorial do vetor 𝒎

𝐷..........................................Matriz de transformação de 𝑆𝑅 para 𝑆𝐶

𝝉𝐶𝑀𝑒𝑥𝑡.......................................Torque externo com relação ao CM

𝒉𝐶𝑀𝐼 .......................................Quantidade de momento angular com relação ao CM escrito na base 𝑆𝐼

𝑘............................................Constante relacionada ao torque interno gerado pelos motores

𝒉𝐶𝑀𝐶 .......................................Quantidade de Momento angular com relação ao CM escrito na base

𝑆𝐶

𝐽𝐶...........................................Matriz de inércia com relação a 𝑆𝐶

𝑚...........................................Massa do quadricóptero

𝑔............................................Aceleração da gravidade

𝜙............................................Roll (ângulo de Euler)

𝜃............................................Pitch (ângulo de Euler)

𝜓...........................................Yaw (ângulo de Euler)

𝒇............................................Força de propulsão gerada pelas hélices

𝒗............................................Velocidade do centro de massa do quadricóptero

�̅�............................................Comando de posição com relação a 𝑆𝐼

�̅�..........................................Comando de atitude segundo a representação por MRP

�̂�............................................Vetor unitário normal ao plano que contém os rotores do

quadricóptero

𝒏𝟏𝟐........................................Projeção de �̂� no plano horizontal

𝑛3..........................................Componente de �̂� em 𝒁𝑹

𝒌𝟏e 𝒌𝟑...................................Constante proporcional do Controlador de Posição

𝒌𝟐e 𝒌𝟒...................................Constante derivativa do Controlador de Posição

𝒌𝒑..........................................Constante proporcional do Controlador de Atitude

𝒌𝒅..........................................Constante derivativa do Controlador de Atitude

(𝑥1, 𝑦1, 𝑧1).............................Coordenadas em 𝑆𝐼 do primeiro ponto no espaço a ser visitado

(𝑥2, 𝑦2, 𝑧2).............................Coordenadas em 𝑆𝐼 do segundo ponto no espaço a ser visitado

(𝑥3, 𝑦3, 𝑧3).............................Coordenadas em 𝑆𝐼 do terceiro ponto no espaço a ser visitado

�̅�(𝑡), �̅�(𝑡), 𝑧̅(𝑡)......................Coordenadas em 𝑆𝐼 do comando temporal de posição

𝑡1............................................Tempo em que o comando de posição chega ao primeiro ponto a ser

visitado

𝑡2............................................Tempo em que o comando de posição chega ao segundo ponto a ser

visitado

𝑡3............................................Tempo em que o comando de posição chega ao terceiro ponto a ser

visitado

𝑣𝑥1..........................................Componente 𝑥 da velocidade de comando no caminho para o

primeiro ponto a ser visitado

𝑣𝑥2..........................................Componente 𝑥 da velocidade de comando no caminho para o segundo

ponto a ser visitado

𝑣𝑥3..........................................Componente 𝑥 da velocidade de comando no caminho para o terceiro

ponto a ser visitado

𝑑1............................................Distância entre o primeiro ponto a ser visitado e a origem

𝑑2............................................Distância entre o primeiro e o segundo ponto a ser visitado

𝑑3............................................Distância entre o segundo e o terceiro ponto a ser visitado

𝑓𝑚í𝑛.........................................Força total mínima gerada pelos rotores

𝑓𝑚í𝑛.........................................Força total máxima gerada pelos rotores

𝛾𝑚á𝑥........................................Máximo ângulo de heading permitido

𝑘𝑎...........................................Constante aerodinâmica referente às pás dos rotores

𝜔𝑚..........................................Velocidade angular nos rotores

Sumário

1 Introdução..........................................................................................................................16

2 Objetivo.............................................................................................................................16

3 Modelagem Dinâmica.......................................................................................................17

3.1 Movimento de Rotação.........................................................................................19

3.2 Movimento de Translação.....................................................................................23

4 Parâmetros Físicos de um Quadricóptero Real.................................................................23

5 Validação das Equações Dinâmicas..................................................................................25

5.1 Funções Utilizadas................................................................................................26

5.2 Testes....................................................................................................................27

6 Sistemas de Controle de Atitude e Posição.......................................................................34

6.1 Controlador de Posição.........................................................................................36

6.2 Controlador de Atitude.........................................................................................40

7 Implementação da Malha Fechada....................................................................................40

8 Resultados e Discussão.....................................................................................................45

9 Conclusões e Sugestões para Trabalhos Posteriores.........................................................49

Referências...........................................................................................................................50

Apêndice – Códigos do MATLAB......................................................................................51

1 Introdução

A motivação deste trabalho provém principalmente da observação das diversas aplicações

que há para quadricópteros. Uma delas é a prospecção de áreas hostis, importante na área militar.

Há também empresas de compras online que utilizam quadricópteros para efetuar a entrega da

mercadoria comprada. Em geral, o transporte autônomo de cargas é corriqueiro, principalmente no

área militar.

Em segundo lugar, há uma motivação proveniente da simples curiosidade de um

Engenheiro que observa um corpo com seis graus de liberdade e deseja compreender sua dinâmica.

O contexto deste trabalho envolve a controlabilidade de um quadricóptero, que possui seis

graus de liberdade. Contudo, há somente quatro variáveis independentes de atuação: a força

exercida por cada uma das quatro hélices. Têm-se, portanto, um sistema subatuado – não é possível

controlar independentemente cada um dos seis graus de liberdade, ou seja, não possível realizar

qualquer manobra com esse corpo.

No entanto, como neste trabalho basta o corpo mover-se de um ponto a outro, é suficiente

poder-se controlar independentemente somente quatro graus de liberdade, que são a posição no

espaço mais o ângulo de heading.

Os controladores de atitude e de posição utilizados são do tipo proporcional-derivativo

(PD) com saturação. Basicamente, o termo proporcional é necessário para tornar o sistema

responsivo, o termo derivativo é importante para suavizar a trajetória, e a saturação impede que o

corpo fique com as hélices invertidas.

O presente texto está organizado da seguinte forma. No Capítulo 3, são definidas as

equações dinâmica; no Capítulo 4, elas serão testadas em ambiente Simulink; no Capítulo 5,

definir-se-á o modelo físico do quadricóptero; no Capítulo 6, serão definidas as equações de

controle; no Capítulo 7, elas serão implementadas em ambiente Simulink; e, no Capítulo 8, os

resultados serão exibidos e discutidos.

2 Objetivo

Por meio deste trabalho, buscou-se compreender toda a modelagem por trás de um

quadricóptero, cuja missão é ir de um ponto a outro autonomamente. Mais especificamente,

pretendeu-se dominar a modelagem física e o controle de posição e atitude do quadricóptero. Por

fim, deseja-se obter uma trajetória suave para uma rota comandada que consiste em segmentos de

reta entre os pontos a serem visitados.

3 Modelagem Dinâmica

Utilizou-se a referência [1] para extrair grande parte das equações dinâmicas.

Definiram-se os sistemas de coordenadas conforme a Fig. 3.1:

Figura 3.1 – Esquema do quadricóptero e definição dos sistemas de coordenas

Sejam: 𝑆𝐼: {𝑋𝐼 , 𝑌𝐼 , 𝑍𝐼}, 𝑆𝑅: {𝑋𝑅 , 𝑌𝑅 , 𝑍𝑅}, 𝑆𝐶 : {𝑋𝐶 , 𝑌𝐶 , 𝑍𝐶}.

Ainda:

𝒓𝑪𝑴 = [

𝑥𝐶𝑀𝑦𝐶𝑀𝑧𝐶𝑀

]

O sistema de coordenadas 𝑆𝑅 é somente uma translação de 𝑆𝐼 e permanece centrado no CM

do quadricóptero.

O sistema de coordenadas 𝑆𝐶 também está centrado no CM, mas ele é fixo no corpo,

rotacionando junto com ele. Mais especificamente, como ilustrado na Fig. 3.1, sua componente 𝑍𝐶

é normal ao plano que contém o quadricóptero, 𝑋𝐶 está no plano que contém o quadricóptero e está

a 45° entre os rotores 1 e 2, e 𝑌𝐶 está no plano que contém o quadricóptero e está a 45° dos rotores

1 e 4.

Note ainda que todos os sistemas de coordenadas são dextrogiros.

A motivação por trás dessa escolha de sistemas de coordenadas é tornar a modelagem física

e as leis de controle o mais simples possível.

Premissas:

1) O veículo é um corpo rígido

A rigor, trata-se de um corpo rígido exceto pelo movimento dos hélices. Não chega a

se tratar de um corpo flexível, mas é mais apropriado dizer que o quadricóptero é um corpo

rígido com movimento interno dos rotores. No entanto, é desprezível a variação da matriz

de inércia devido ao movimento das hélices.

2) 𝑺𝑰 é inercial, pois o movimento da terra é desprezível com relação ao do

quadricóptero

É plausível desprezar o movimento da Terra porque ele é muito maior que a

distância percorrida pelo quadricóptero, que neste trabalho é menor que 40𝑚, porque a

altitude do movimento é muito baixa e porque a velocidade do corpo é pequena quando

comparada a da Terra.

3) O veículo tem distribuição de massa simétrica com relação a 𝑺𝑪

Esta premissa é importante porque simplifica as equações dinâmicas de rotação. O

prejuízo é mínimo uma vez que os termos que não são da diagonal principal da matriz de

inércia são menores que 0,4% dos termos da diagonal principal, como será mostrado no

capítulo 4.

A premissa 3 é equivalente a dizer que 𝑆𝐶 define os eixos principais de inércia do

quadricóptero.

3.1 Movimento de Rotação

A atitude de 𝑆𝐶 com respeito a 𝑆𝑅 será representada pelo vetor dos Parâmetros de

Rodrigues Modificados (MRP):

𝒎 = [

𝑚1

𝑚2

𝑚3

]

Por praticidade, utilizaram-se as expressões “atitude do corpo” ou “atitude do

quadricóptero” para se referir à rotação de 𝑆𝐶 com respeito a 𝑆𝑅.

A velocidade angular, escrita na base 𝑆𝐶, de 𝑆𝐶 com relação a 𝑆𝑅 será denotada por:

𝝎 = [

𝜔𝑥𝜔𝑦𝜔𝑧]

Significado Geométrico dos MRP

Bastam 3 parâmetros independentes para se definir a rotação de um sistema de coordenadas

com relação a outro, ou seja, atitude.

Uma dessas formas é usar a representação de Eixo e Ângulo de Euler, que verificou que a

atitude de qualquer sistema de coordenas é resultado da rotação de um ângulo 𝜑 de outro sistema

de coordenas com relação a um vetor unitário �̂�. De fato, têm-se aqui somente três parâmetros

independentes: o ângulo𝜑 e duas componentes de �̂�, já que ele é um vetor unitário.

Figura 3.2 – Rotação do sistema de coordenadas segundo Eixo e Ângulo de Euler

A representação por MRP, que também possui três parâmetros independentes, é outra

forma mais concisa de se trabalhar, já que facilita a escrita das equações cinemáticas, como será

visto adiante. A partir da representação por Eixo e Ângulo de Euler, é possível obter os MRPs a

partir da seguinte relação:

𝒎 = 𝑡𝑔 (𝜑

4) �̂� (𝟑. 𝟏)

Haveria uma singularidade se 𝜑 = 360°, por exemplo. Contudo, isso não ocorreria em

nenhum momento da operação do quadricóptero, cuja única missão é ir de um ponto a outro. Se o

quadricóptero precisa fazer um movimento lateral, 𝜑 ficaria no intervalo ] − 90°, 90°[, se ele

precisa subir ou descer, nem precisa girar. Por outro lado, uma manobra em que o quadricóptero

gira em torno de si mesmo indefinidamente não poderia ser realizada, mas ela não está no escopo

deste trabalho.

Por fim, vale ressaltar que os Ângulos de Euler não foram usados nas equações dinâmicas

porque, neste caso, haveria singularidades em 90° e o esforço computacional seria maior por haver

funções trigonométricas.

Equação Cinemática

Pode-se mostrar que a evolução temporal da atitude do quadricóptero é regida pela seguinte

EDO:

�̇� =1

4{(1 −𝒎𝑡𝒎)𝐼3 + 2[𝑚×] + 2𝒎𝒎

𝑡}𝝎 (𝟑. 𝟐)

Onde: [𝑚×] = [

0 −𝑚3 𝑚2

𝑚3 0 −𝑚1

−𝑚2 𝑚1 0].

A matriz de transformação de 𝑆𝑅 para 𝑆𝐶 é dada por:

𝐷 = 𝐼3 +1

(1 +𝒎𝑡𝒎)2{8[𝑚×]

2 − 4(1 −𝒎𝑡𝒎)[𝑚×]} (𝟑. 𝟑)

Apesar das 𝐸𝑞. (3.2) e 𝐸𝑞. (3.3) apresentarem uma forma relativamente simples, quando

comparadas às equações equivalentes para Ângulos de Euler, por exemplo, a demonstração delas

requer grande esforço algébrico.

Agora fica clara mais uma vantagem de se usar os MRPs: não há funções trigonométricas,

o que deixa a simulação mais rápida e evita o aparecimento de singularidades. É possível verificar

que, de fato, não há singularidades na 𝐸𝑞. (3.3):

𝒎𝑡𝒎 = 𝑚12 +𝑚2

2 +𝑚32 = [𝑡𝑔 (

𝜑

4) 𝑒1]

2

+ [𝑡𝑔 (𝜑

4) 𝑒2]

2

+ [𝑡𝑔 (𝜑

4) 𝑒3]

2

= 𝑡𝑔2 (𝜑

4) (𝑒1

2 + 𝑒22 + 𝑒3

2) = 𝑡𝑔2 (𝜃

4) ≥ 0 ⇒ (1 +𝒎𝑡𝒎)2 ≥ 1

Portanto, o denominador da 𝐸𝑞. (3.3) nunca se anula.

Equação Dinâmica

Para escrever algumas equações dinâmicas de rotação, foram observadas notações da

referência [4]. Adota-se aqui a formulação Newton-Euler.

Pela 2a Lei de Euler, a equação dinâmica de rotação de um corpo rígido é

∑𝝉𝐶𝑀𝑒𝑥𝑡 =

𝑑𝒉𝐶𝑀𝐼

𝑑𝑡 (𝟑. 𝟒)

Observando-se a Fig. 3.1 e notando-se que 𝒇𝒊 é perpendicular ao plano que contém o

quadricóptero, fica fácil calcular o lado esquerdo da 𝐸𝑞. (3.4):

∑𝝉𝐶𝑀𝑒𝑥𝑡 = [

(+𝑓1 − 𝑓2 − 𝑓3 + 𝑓4)ℓ√2 2⁄

(−𝑓1 − 𝑓2 + 𝑓3 + 𝑓4)ℓ√2 2⁄

(+𝑓1 − 𝑓2 + 𝑓3 − 𝑓4)𝑘

]

𝐶

= [

𝜏𝑥𝜏𝑦𝜏𝑧]

𝐶

(𝟑. 𝟓)

Os torques da 𝐸𝑞. (3.5) foram definidos positivos segundo a regra-da-mão-direita.

A constante 𝑘 possui dimensão de distância. Segundo o princípio da conservação da

quantidade de momento angular, o torque na direção 𝑍𝐶 aparece para compensar os desequilíbrio

de torques internos nos rotores. Agora é possível entender que as pás na Fig. 3.1 não giram todas

no mesmo sentido porque, se não, haveria um yaw indesejável. Assim, as hélices dos rotores 1 e 3

giram no sentido horário, enquanto as dos rotores 2 e 4 giram no sentido anti-horário.

Para que se possa escrever 𝒉𝐶𝑀𝐼 na base 𝑆𝐼, é preciso o seguinte ajuste por 𝑆𝐶 não ser um

sistema de coordenadas inercial:

𝑑𝒉𝐶𝑀𝐼

𝑑𝑡=𝑑𝒉𝐶𝑀

𝐶

𝑑𝑡+ 𝝎 × 𝒉𝐶𝑀

𝐶 (𝟑. 𝟔)

Pode-se escrever 𝒉𝐶𝑀𝐶 como:

𝒉𝐶𝑀𝐶 = 𝐽𝐶𝝎 (𝟑. 𝟕)

Onde: 𝐽𝐶 = [

𝐽𝑥𝑥 −𝐽𝑥𝑦 −𝐽𝑥𝑧−𝐽𝑦𝑥 𝐽𝑦𝑦 −𝐽𝑦𝑧−𝐽𝑧𝑥 −𝐽𝑧𝑦 𝐽𝑧𝑧

] é obtida com respeito a 𝑆𝐶.

Substituindo a 𝐸𝑞. (3.7) na 𝐸𝑞. (3.6), vem:

𝑑𝒉𝐶𝑀𝐼

𝑑𝑡= 𝐽𝐶�̇� + 𝝎 × (𝐽𝐶𝝎) (𝟑. 𝟖)

Substituindo as 𝐸𝑞. (3.8) e 𝐸𝑞. (3.5) na 𝐸𝑞. (3.4), e considerando 𝐽𝐶 = 𝑑𝑖𝑎𝑔(𝐽𝑥, 𝐽𝑦, 𝐽𝑧),

vem:

{

�̇�𝑥 =

𝐽𝑦 − 𝐽𝑧

𝐽𝑥𝜔𝑦𝜔𝑧 +

1

𝐽𝑥𝜏𝑥

�̇�𝑦 =𝐽𝑧 − 𝐽𝑥𝐽𝑦

𝜔𝑥𝜔𝑧 +1

𝐽𝑦𝜏𝑦

�̇�𝑧 =𝐽𝑥 − 𝐽𝑦

𝐽𝑧𝜔𝑥𝜔𝑦 +

1

𝐽𝑧𝜏𝑧

(𝟑. 𝟗)

3.2 Movimento de Translação

Pela 2a Lei de Newton:

𝑭𝑅 = 𝑚𝒂𝐶𝑀 (𝟑. 𝟏𝟎)

Escrevendo a 𝐸𝑞. (3.10) no referencial 𝑆𝐼:

[00

−𝑚𝑔] + 𝐷𝑡 [

00

𝑓1 + 𝑓2 + 𝑓3 + 𝑓4

] = 𝑚 [

�̈�𝐶𝑀�̈�𝐶𝑀�̈�𝐶𝑀

] (𝟑. 𝟏𝟏)

Por fim, a dinâmica do quadricóptero é definida pelos sistemas de EDOs das 𝐸𝑞. (3.2),

𝐸𝑞. (3.9) e 𝐸𝑞. (3.11). As incógnitas são a posição do CM do corpo, sua velocidade angular e sua

atitude.

4 Parâmetros Físicos de um Quadricóptero Real

A simulação será feita usando dados do quadricóptero Gyro200 ED 4X da Gyrofly [2].

Esse equipamento está ilustrado na Fig. 4.1. A Tabela 4.1 lista as estimativas dos parâmetros do

quadricóptero em questão fornecidas pelo fabricante.

Figura 4.1 – Desenho do quadricóptero em CAD

Massa 𝒎 (𝒌𝒈) Braço 𝒍 (𝒎)

1,2 0,277

Tabela 4.1 – Parâmetros do Quadricóptero

Extraíram-se também os momentos de inércia em 𝑔 ∙ 𝑚𝑚2, medidos com relação ao

sistema de coordenadas fixo no corpo:

Tabela 4.2 – Momentos de inércia do corpo em 𝒈.𝒎𝒎𝟐

É razoável assumir que a matriz de inércia da Tabela 4.2 é diagonal porque os valores fora

da diagonal principal são menores que 0,4% dos valores da diagonal principal. Isto comprova a

premissa 3 do Capítulo 3 e permite que os valores da diagonal principal da Tabela 4.1 sejam

aproveitados neste trabalho:

𝐽𝑥 = 𝐼𝑥𝑥 = 0,02259538940 𝑘𝑔 ∙ 𝑚2

𝐽𝑦 = 𝐼𝑦𝑦 = 0,02218007586 𝑘𝑔 ∙ 𝑚2

𝐽𝑧 = 𝐼𝑧𝑧 = 0,03761572017 𝑘𝑔 ∙ 𝑚2

Por mais que, em [2], o eixo 𝑍𝐶 aponte para baixo, enquanto neste trabalho ele aponta para

cima, conforme a Fig. 3.1, o termo 𝐼𝑧𝑧 =∭ 𝑧2𝑑𝑚𝑉

não mudaria, pois o termo 𝑧2 anula o efeito

do sinal da componente 𝑧 dos elementos de massa.

5 Validação das Equações Dinâmicas

Como próximo passo, poder-se-ia implementar a malha de controle e depois realizar testes.

O problema é que, caso houvesse algum problema durante a simulação da malha fechada, seria

muito difícil identificar se o problema está na dinâmica do quadricóptero ou na malha de controle.

Dessa forma, optou-se por testar a dinâmica do corpo em primeiro lugar. Só quando se

tivesse certeza que ela está funcionando corretamente, implementar-se-ia a malha de controle.

Neste trabalho, foi adotado 𝑔 = 9,81 𝑚/𝑠2. Na sequência, veem-se os blocos em ambiente

Simulink em que foram realizados os testes.

Figura 5.1 – Ambiente Simulink para testes da dinâmica

No apêndice encontram-se todos os códigos utilizados em MATLAB.

5.1 Funções Utilizadas

A seguir, serão explicados os principais blocos da Fig. 5.1.

dinamica

Foram implementadas as 𝐸𝑞. (3.2), 𝐸𝑞. (3.9) e 𝐸𝑞. (3.11) como equações dinâmicas e a

𝐸𝑞. (3.3) como equação estática. O arquivo é do tipo S-Function e possui um formato característico

para resolver sistemas de EDOs de primeiro grau, como é o caso deste trabalho.

m2angulosEuler

Apesar de os MRPs serem a forma mais simples de representar atitude nas equações

dinâmicas, essa representação não é a melhor maneira para visualizar a atitude do corpo. Dessa

forma, primeiramente converteram-se os MRP para a matriz de atitude através da 𝐸𝑞. (3.3), e, em

seguida, converteu-se a matriz de atitude 𝐷 para a representação por Ângulos de Euler.

Os Ângulos de Euler são: roll 𝜙, pitch 𝜃 e yaw 𝜓. A ordem utilizada foi a seguinte:

primeiro há um roll sobre o eixo 𝑋𝐶, depois, a partir da configuração obtida, um pitch sobre o

eixo𝑌𝐶 e, por último, a partir da configuração obtida, um yaw sobre o eixo 𝑍𝐶 . Essa ordem é

caracterizada como 1-2-3.

As fórmulas de conversão de 𝐷 para Ângulos de Euler são:

𝜃 = −sen−1(−𝐷31) (𝟓. 𝟏)

𝜙 = −𝑡𝑔−1(𝐷32, 𝐷33) (𝟓. 𝟐)

𝜓 = −𝑡𝑔−1(𝐷21, 𝐷11) (𝟓. 𝟑)

forces

Esse é o bloco que impõe as forças nos rotores para que determinado teste seja realizado.

5.2 Testes

Para testar por completo se a dinâmica foi implementada corretamente, foram realizados

quatro testes: pairar, subir, movimento lateral e yaw puro.

Pairar

𝑓1 = 𝑓2 = 𝑓3 = 𝑓4 =𝑃

4 (𝟓. 𝟒)

Condições iniciais: 𝒓𝑪𝑴 = 0 𝑚, 𝒗𝑪𝑴 = 0𝑚 𝑠⁄ ,𝒎 = 0,𝝎 = 0 𝑟𝑎𝑑 𝑠⁄

Resultados:

Fig. 5.2 – Pairar | 𝒙𝑪𝑴 Fig. 5.3 – Pairar | 𝒚𝑪𝑴 Fig. 5.4 – Pairar | 𝒛𝑪𝑴

Fig. 5.5 – Pairar | �̇�𝑪𝑴 Fig. 5.6 – Pairar | �̇�𝑪𝑴 Fig. 5.7 – Pairar | �̇�𝑪𝑴

Fig. 5.8 – Pairar | 𝜙 Fig. 5.9 – Pairar | 𝜽 Fig. 5.10 – Pairar | 𝝍

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

xC

M (

m)

t (s)

xCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10y

CM

(m

)

t (s)

yCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

z CM

(m

)

t (s)

zCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vxC

M (

m/s

)

t (s)

vxCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vyC

M (

m/s

)

t (s)

vyCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vzC

M (

m/s

)

t (s)

vzCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

Ro

ll (ᵒ)

t (s)

Roll

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

Pit

ch (ᵒ)

t (s)

Pitch

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

Ya

w (ᵒ)

t (s)

Yaw

Fig. 5.11 – Pairar | 𝝎𝒙 Fig. 5.12 – Pairar | 𝝎𝒚 Fig. 5.13 – Pairar | 𝝎𝒛

Subir

𝑓1 = 𝑓2 = 𝑓3 = 𝑓4 =𝑃

2 (𝟓. 𝟓)

Há uma força resultante vertical igual a 𝑃 que produzirá uma aceleração constante,

velocidade linear e posição parabólica, como mostrado nas Fig. 5.14 a 5.19 abaixo.

Condições iniciais: 𝒓𝑪𝑴 = 0 𝑚, 𝒗𝑪𝑴 = 0𝑚 𝑠⁄ ,𝒎 = 0,𝝎 = 0 𝑟𝑎𝑑 𝑠⁄

Resultados:

Fig. 5.14 – Subir | 𝒙𝑪𝑴 Fig. 5.15 – Subir | 𝒚𝑪𝑴 Fig. 5.16 – Subir | 𝒛𝑪𝑴

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωx (

rad

/s)

t (s)

ωx

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωy (

rad

/s)

t (s)

ωy

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωz

(ra

d/s

)

t (s)

ωz

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

xC

M (

m)

t (s)

xCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

yC

M (

m)

t (s)

yCM

0

100

200

300

400

500

600

0 1 2 3 4 5 6 7 8 9 10

z CM

(m

)

t (s)

zCM

Fig. 5.17 – Subir | �̇�𝑪𝑴 Fig. 5.18 – Subir | �̇�𝑪𝑴 Fig. 5.19 – Subir | �̇�𝑪𝑴

Fig. 5.20 – Subir | 𝜙 Fig. 5.21 – Subir | 𝜽 Fig. 5.22 – Subir | 𝝍

Fig. 5.23 – Subir | 𝝎𝒙 Fig. 5.24 – Subir | 𝝎𝒚 Fig. 5.25 – Subir | 𝝎𝒛

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vxC

M (

m/s

)

t (s)

vxCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vyC

M (

m/s

) t (s)

vyCM

0

20

40

60

80

100

120

0 1 2 3 4 5 6 7 8 9 10

vzC

M (

m/s

)

t (s)

vzCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

Ro

ll (ᵒ)

t (s)

Roll

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

Pit

ch (ᵒ)

t (s)

Pitch

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

Ya

w (ᵒ)

t (s)

Yaw

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωx (

rad

/s)

t (s)

ωx

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωy (

rad

/s)

t (s)

ωy

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωz

(ra

d/s

)

t (s)

ωz

Lateral

Desejou-se inclinar o corpo de 30° na direção, posicionando-se o rotor 1 à frente. Assim, o

eixo de rotação deve estar na direção [

𝑋𝐶𝑌𝐶𝑍𝐶

] = [−1+10]. Assim, pela 𝐸𝑞. (3.1), vem:

𝒎(0) = 𝑡𝑔 (30°

4) [−√2 2⁄

+√2 2⁄0

] (𝟓. 𝟔)

As forças são tais que há uma resultante horizontal de módulo 𝑃√3

3.

𝑓1 = 𝑓2 = 𝑓3 = 𝑓4 =𝑃√3

6 (𝟓. 𝟕)

Como há uma resultante horizontal constante, velocidade é linear e a posição é parabólica,

como se pode verificar nas Fig. 5.26 a 5.31 abaixo.

Condições iniciais: 𝒓𝑪𝑴 = 0 𝑚, 𝒗𝑪𝑴 = 0𝑚 𝑠⁄ ,𝒎 = [−0,0931+0,0931

0] ,𝝎 = 0 𝑟𝑎𝑑 𝑠⁄

Resultados:

Fig. 5.26 – Lateral | 𝒙𝑪𝑴 Fig. 5.27 – Lateral | 𝒚𝑪𝑴 Fig. 5.28 – Lateral | 𝒛𝑪𝑴

0

50

100

150

200

250

0 1 2 3 4 5 6 7 8 9 10

xC

M (

m)

t (s)

xCM

0

50

100

150

200

250

0 1 2 3 4 5 6 7 8 9 10

yC

M (

m)

t (s)

yCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

z CM

(m

)

t (s)

zCM

Fig. 5.29 – Lateral | �̇�𝑪𝑴 Fig. 5.30 – Lateral | �̇�𝑪𝑴 Fig. 5.31 – Lateral | �̇�𝑪𝑴

Fig. 5.32 – Lateral | 𝜙 Fig. 5.33 – Lateral | 𝜽 Fig. 5.34 – Lateral | 𝝍

Fig. 5.35 – Lateral | 𝝎𝒙 Fig. 5.36 – Lateral | 𝝎𝒚 Fig. 5.37 – Lateral | 𝝎𝒛

0

10

20

30

40

50

0 1 2 3 4 5 6 7 8 9 10

vxC

M (

m/s

)

t (s)

vxCM

0

10

20

30

40

50

0 1 2 3 4 5 6 7 8 9 10

vyC

M (

m/s

) t (s)

vyCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vzC

M (

m/s

)

t (s)

vzCM

-25

-20

-15

-10

-5

0

0 1 2 3 4 5 6 7 8 9 10

Ro

ll (ᵒ)

t (s)

Roll

0

5

10

15

20

25

0 1 2 3 4 5 6 7 8 9 10

Pit

ch (ᵒ)

t (s)

Pitch

0

1

2

3

4

5

0 1 2 3 4 5 6 7 8 9 10

Ya

w (ᵒ)

t (s)

Yaw

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωx (

rad

/s)

t (s)

ωx

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωy (

rad

/s)

t (s)

ωy

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

ωz

(ra

d/s

)

t (s)

ωz

Rotação pura

𝑓1 = 𝑓3 =3𝑃

8, 𝑓2 = 𝑓4 =

𝑃

8 (𝟓. 𝟖)

A escolha das forças na 𝐸𝑞. (5.8) garante que (i) a resultante das forças seja nula e (ii) haja

um desequilíbrio de quantidade de momento angular que faz o corpo girar no sentido positivo. Tais

eventos podem ser verificados nas Fig. 5.38 a 5.49 abaixo.

Condições iniciais: 𝒓𝑪𝑴 = 0 𝑚, 𝒗𝑪𝑴 = 0𝑚 𝑠⁄ ,𝒎 = 0,𝝎 = 0 𝑟𝑎𝑑 𝑠⁄

Resultados:

Fig. 5.38 – Yaw | 𝒙𝑪𝑴 Fig. 5.39 – Yaw | 𝒚𝑪𝑴 Fig. 5.40 – Yaw | 𝒛𝑪𝑴

Fig. 5.41 – Yaw | �̇�𝑪𝑴 Fig. 5.42 – Yaw | �̇�𝑪𝑴 Fig. 5.43 – Yaw | �̇�𝑪𝑴

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

xC

M (

m)

t (s)

xCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

yC

M (

m)

t (s)

yCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10z C

M (

m)

t (s)

zCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vxC

M (

m/s

)

t (s)

vxCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vyC

M (

m/s

)

t (s)

vyCM

0,0

0,2

0,4

0,6

0,8

1,0

0 1 2 3 4 5 6 7 8 9 10

vzC

M (

m/s

)

t (s)

vzCM

Fig. 5.44 – Yaw | 𝜙 Fig. 5.45 – Yaw | 𝜽 Fig. 5.46 – Yaw | 𝝍

Fig. 5.47 – Yaw | 𝝎𝒙 Fig. 5.48 – Yaw | 𝝎𝒚 Fig. 5.49 – Yaw | 𝝎𝒛

Ocorre um erro na simulação em 𝑡 = 3,1 𝑠, pois, neste instante, o corpo deu uma volta

completa em torno de 𝑍𝐶 . Essa singularidade estava prevista, conforme discutido no Capítulo 3.1.

Como todos os resultados mostraram a saída esperada, pôde-se concluir que a dinâmica foi

implementada corretamente. Com isso, prosseguiu-se para a questão do controle.

6 Sistemas de Controle de Atitude e Posição

É importante, antes de mais nada, entender a importância e a necessidade de se utilizar um

sistema de controle.

Primeiro, o quadricóptero está sujeito a distúrbios inesperados que podem tirá-lo de uma

rota predeterminada. O exemplo mais claro disso são os ventos, totalmente aleatórios. Podem-se ter

-1,0

-0,5

0,0

0,5

1,0

0 1 2 3 4 5 6 7 8 9 10Ro

ll (ᵒ)

t (s)

Roll

-1,0

-0,5

0,0

0,5

1,0

0 1 2 3 4 5 6 7 8 9 10Pit

ch (ᵒ)

t (s)

Pitch

-200

-100

0

100

200

0 1 2 3 4 5 6 7 8 9 10Ya

w (ᵒ)

t (s)

Yaw

-1,0

-0,5

0,0

0,5

1,0

0 1 2 3 4 5 6 7 8 9 10

ωx (

rad

/s)

t (s)

ωx

-1,0

-0,5

0,0

0,5

1,0

0 1 2 3 4 5 6 7 8 9 10

ωy (

rad

/s)

t (s)

ωy

0

1

2

3

4

5

0 1 2 3 4 5 6 7 8 9 10

ωz

(ra

d/s

)

t (s)

ωz

ainda atritos e vibrações mecânicas que não foram modeladas. Assim, é necessário um sistema de

controle que procure manter o quadricóptero sempre no ponto de operação desejado e que seja

robusto o bastante para não deixar o quadricóptero se desestabilizar e cair.

Segundo, as próprias equações e/ou parâmetros do modelo podem não representar o

sistema físico de forma exata. Muitas vezes as EDOs ou EDPs que modelam um sistema são

linearizadas (não é o caso deste trabalho), surgindo a necessidade de “se proteger” disso através de

uma malha de controle. É possível também que as medidas do quadricóptero não estejam 100%

corretas: a massa pode ter sido medida numa balança descalibrada, as medidas do corpo podem não

ter sido aferidas com tanta precisão, etc. Todos esses erros podem também se acumular na matriz

de inércia.

Vê-se, portanto, que a malha de controle é essencial para fazer um sistema funcionar com

tantas imperfeições e imprevisibilidades. A seguir, vê-se um esquema simplificado das malhas

fechadas:

Figura 6.1 – Esquema simplificado das malhas fechadas do controlador do quadricóptero

A malha da Fig. 6.1 opera da seguinte maneira: o Controlador de Posição compara a

posição atual com a comandada, calculando um erro. Com base neste erro, determina a intensidade

do thrust total e a atitude necessária para atingir a posição desejada. Mas ainda falta efetivamente

levar o quadricóptero da atitude atual à comandada, e quem faz isso é o Controlador de Atitude.

Similarmente ao Controlador de Posição, ele calcula a diferença entre a atitude atual e a desejada,

e, com base neste erro, envia comandos de torque ao quadricóptero. O thrust e torque determinados

pelos controladores causam uma modificação no estado do corpo. O novo estado, agora mais

próximo da posição desejada, é então realimentado aos controladores para se repetir o ciclo.

Inicialmente, o alto valor do erro produzirá comandos de thrust e torque altos, pois se está

longe da posição desejada. Conforme o quadricóptero se aproxima do ponto desejado, o erro

diminui e, com ele, os comandos de força nos rotores. Quando se está muito próximo do ponto

desejado, após algumas iterações do ciclo, o erro é praticamente zero, assim como as forças de

atuação, que não mais alteram o estado do quadricóptero. Neste momento, o ciclo para de iterar.

Há duas malhas de controle: uma para posição e outra para atitude, que depende da saída

do controlador de posição. Essa relação é intuitiva: qualquer mudança horizontal na posição requer

uma determinada inclinação do quadricóptero. Em outras palavras, o Controlador de Posição

determina a atitude necessária, pois o movimento de translação horizontal depende da atitude. O

Controlador de Atitude, por sua vez, é responsável por fazer a atitude atual igualar-se à atitude

necessária para se chegar a determinado ponto. Fica clara, portanto, a necessidade de haver duas

malhas de controle.

6.1 Controlador de Posição

As leis do Controlador de Posição foram extraídas da referência [3].

No esquema a seguir, �̂� = 𝒁𝑪/|𝒁𝑪| e 𝒏𝟏𝟐 e 𝒏𝟑 são suas projeções no plano horizontal e em

𝒁𝑹, respectivamente.

Figura 6.2 – Esquema auxiliar para determinação das leis de controle

Pela Figura 6.2, é possível visualizar como ocorre o movimento puramente lateral: inclina-

se o quadricóptero na direção que se deseja ir através de um aumento do thrust proveniente da parte

posterior do quadricóptero. Com isso, é necessário aumentar o thrust gerado pelos rotores, de

forma a equilibrar o peso e manter o movimento na horizontal. O movimento puramente vertical é

bem simples: basta aumentar ou diminuir o thrust.

Quando o movimento é horizontal e vertical, há uma combinação dos mecanismos

descritos acima. Pela diferença dos mecanismos, fica claro que é preciso separar o controle em

duas partes: movimento vertical e horizontal. Será comprovado a seguir que o thrust 𝑓 é

responsável por comandar o movimento vertical, enquanto a componente 𝑛12 de 𝒏 é responsável

por comandar o movimento horizontal.

Movimento Vertical

Para se determinar a lei de controle do movimento vertical, analisou-se a equação do

movimento vertical de translação. Assim, escreveu-se a terceira componente da 𝐸𝑞. (3.11):

−𝑔 + 𝑓𝑛3𝑚= �̈�3 (𝟔. 𝟏)

Com:

𝑓 = 𝑓1 + 𝑓2 + 𝑓3 + 𝑓4 (𝟔. 𝟐)

Olhando para a 𝐸𝑞. (6.1), vê-se que se pode controlar �̈�3através de 𝑛3 = 𝐷33, mas isso

alteraria o movimento horizontal, conforme discutido acima. Assim, controla-se �̈�3através de 𝑓.

Vale lembrar que a matriz de atitude 𝐷 é obtida através da 𝐸𝑞. (3.3).

A escolha de 𝑓 deve ser tal que a equação resultante seja da forma:

�̈�3 + 𝑘2�̇�3 + 𝑘1𝑟3 = 𝑘1�̅�3 (𝟔. 𝟑)

É fácil ver que 𝑟3 = �̅�3 é solução da 𝐸𝑞. (6.3), que é justamente o que se deseja obter.

Trata-se de um controle PD, onde 𝑘1 e 𝑘2 são constantes a serem determinadas. Através das

𝐸𝑞. (6.1) e 𝐸𝑞. (6.3), vê-se que 𝑓 deve obedecer à seguinte lei:

𝑓 =𝑚

𝑛3[𝑔 − 𝑘1(𝑟3 − �̅�3) − 𝑘2�̇�3] (𝟔. 𝟒)

Não é preciso preocupar-se com a singularidade 𝑛3 = 0 porque, conforme será visto a

seguir, a saturação definida para o movimento horizontal garante que a projeção de 𝒏 em 𝒁𝑹 nunca

será nula.

A 𝐸𝑞. (6.4) não impõe limites ao thrust total, enquanto, na prática, esses limites existem. É

preciso impor um 𝑓𝑚í𝑛 porque, caso 𝑓 ≅ 0 durante voo, seria muito dispendioso acionar novamente

o sistema eletrônico, o que pode dar margem a acidentes. Para o quadricóptero descer, por

exemplo, teoricamente poder-se-ia ter 𝑓 = 0, mas na prática isso não seria uma boa estratégia. É

preciso também haver um limite superior, 𝑓𝑚á𝑥, porque os motores tem um limite de rotação que

podem impor. Portanto, a 𝐸𝑞. (6.4) deve ter a seguinte forma:

𝑓 = {

𝑓𝑚í𝑛, 𝑠𝑒 𝛼3 < 𝑓𝑚í𝑛𝛼3, 𝑠𝑒 𝑓𝑚í𝑛 ≤ 𝛼3 ≤ 𝑓𝑚á𝑥𝑓𝑚á𝑥, 𝑠𝑒 𝛼3 > 𝑓𝑚á𝑥

(𝟔. 𝟓)

Com:

𝛼3 =𝑚

𝑛3[𝑔 − 𝑘1(𝑟3 − �̅�3) − 𝑘2�̇�3] (𝟔. 𝟔)

Movimento Horizontal

Para se determinar a lei de controle do movimento horizontal, analisou-se a equação do

movimento horizontal translação. Assim, escreveu-se a primeira e segunda componente da

𝐸𝑞. (3.11):

𝑓�̅�𝟏𝟐𝑚

= �̈�𝟏𝟐 (𝟔. 𝟕)

Não é possível controlar 𝒓𝟏𝟐 através de 𝑓, pois 𝑓 não altera a direção de 𝒓𝟏𝟐. Assim,

controla-se 𝒓𝟏𝟐 através de 𝒏𝟏𝟐.

A escolha de 𝑓 deve ser tal que a equação resultante seja da forma:

�̈�𝟏𝟐 + 𝑘4�̇�𝟏𝟐 + 𝑘3𝒓𝟏𝟐 = 𝑘3�̅�𝟏𝟐 (𝟔. 𝟖)

É fácil ver que 𝒓𝟏𝟐 = �̅�𝟏𝟐 é solução da 𝐸𝑞. (6.8), que é justamente o que se deseja obter.

Trata-se de um controle PD, onde 𝑘3 e 𝑘4 são constantes a serem determinadas. Através das

𝐸𝑞. (6.7) e 𝐸𝑞. (6.8), vê-se que �̅�𝟏𝟐 deve obedecer à seguinte lei:

�̅�𝟏𝟐 = −𝑚

𝑓[𝑘3(𝒓𝟏𝟐 − �̅�𝟏𝟐) + 𝑘4�̇�𝟏𝟐] (𝟔. 𝟗)

Não é preciso preocupar-se com a singularidade 𝑓 = 0 porque, conforme foi visto, a

saturação definida para o movimento vertical garante que 𝑓 ≥ 𝑓𝑚í𝑛 > 0.

A 𝐸𝑞. (6.9) não impõe limites a �̅�𝟏𝟐, enquanto, na prática, esses limites devem existir. É

preciso limitar |�̅�𝟏𝟐| = 𝑠𝑒𝑛(𝛾𝑚á𝑥), conforme a Fig. 6.2. Caso contrário, 𝛾𝑚á𝑥 se aproximaria de

90°. Nessa situação, o plano que contém os rotores estaria na vertical, e o corpo cairia, pois a única

força vertical seria o peso. Portanto, a 𝐸𝑞. (6.9) deve ter a seguinte forma:

�̅�𝟏𝟐 = {

𝜶𝟏𝟐, 𝑠𝑒 |𝜶𝟏𝟐| ≤ 𝑠𝑒𝑛(𝛾𝑚á𝑥)

𝑠𝑒𝑛(𝛾𝑚á𝑥)𝜶𝟏𝟐|𝜶𝟏𝟐|

, 𝑠𝑒 |𝜶𝟏𝟐| > 𝑠𝑒𝑛(𝛾𝑚á𝑥) (𝟔. 𝟏𝟎)

Com:

𝜶𝟏𝟐 = −𝑚

𝑓[𝑘3(𝒓𝟏𝟐 − �̅�𝟏𝟐) + 𝑘4�̇�𝟏𝟐] (𝟔. 𝟏𝟏)

Contudo, como se vê na Fig. 6.1, a saída do Controlador de Posição é �̅�, e não �̅�𝟏𝟐.

Fazendo referência à Fig. 3.2, para obter �̅� utiliza-se a 𝐸𝑞. (3.1):

�̅� = 𝑡𝑔 (�̅�

4) �̂� (𝟔. 𝟏𝟐)

Onde:

�̂� =𝒁𝑹 × �̅�𝟏𝟐|𝒁𝑹 × �̅�𝟏𝟐|

(𝟔. 𝟏𝟑)

𝑠𝑒𝑛(�̅�) = |�̅�𝟏𝟐| (𝟔. 𝟏𝟒)

6.2 Controlador de Atitude

A lei de controle para o Controlador de Atitude é relativamente simples:

𝝉 = 𝑘𝑝(�̅� −𝒎) − 𝑘𝑑𝝎 (𝟔. 𝟏𝟓)

7 Implementação da Malha Fechada

Abaixo vêem-se os blocos utilizados no ambiente Simulink e um esquema para a malha

fechada:

Figura 7.1 – Malha fechada do sistema implementada no Simulink

Figura 7.2 – Esquema completo das malhas fechadas do controlador do quadricóptero

Todos os respectivos códigos encontram-se no Apêndice no final deste trabalho. A seguir

serão explicados os blocos da Fig. 7.1 que ainda não foram mencionados.

Em um primeiro momento, era utilizada como entrada uma função degrau ao invés de uma

função rampa. O problema da função degrau é que, quando o comando de posição está longe da

posição atual, o quadricóptero é submetido a uma alta velocidade, podendo sofrer um overshoot.

Além disso, não há neste trabalho a necessidade do quadricóptero mover-se rapidamente.

No lugar de um degrau, foi utilizada, portanto, a função rampas. Dessa forma, há uma

aproximação gradual da posição desejada. Ainda, com a rampa, é possível definir a velocidade

aproximada do corpo – é simplesmente a inclinação da rampa. Para a simulação final, o

quadricóptero deveria sair do ponto (0,0,0), ir para (𝑥1, 𝑦1, 𝑧1) = (7,8,9), depois para

(𝑥2, 𝑦2, 𝑧2) = (10,4,5) e retornar para o ponto de início (𝑥3, 𝑦3, 𝑧3) = (0,0,0), fazendo todo o

percurso a uma velocidade próxima de 𝑣 = 1 𝑚/𝑠, comum entre aplicações de quadricópteros.

Foi criada a função rampasTeste para verificar se a função rampas insere os comandos de

posição corretos. As saídas da função rampasTeste estão ilustradas nas Fig. 7.3 à Fig. 7.5:

Fig. 7.3 – Comando de 𝒙𝑪𝑴 Fig. 7.4 – Comando de 𝒚𝑪𝑴 Fig. 7.5 – Comando de 𝒛𝑪𝑴

Abaixo consta a equação utilizada para se definir os segmentos de rampa �̅� = �̅�(𝑡), sendo

os outros segmentos de rampa �̅� = �̅�(𝑡) e 𝑧̅ = 𝑧̅(𝑡) análogos.

�̅�(𝑡) =

{

𝑣𝑥1𝑡, 𝑡 < 𝑡1

𝑥1 + 𝑣𝑥2(𝑡 − 𝑡1), 𝑡1 ≤ 𝑡 < 𝑡2𝑥2 + 𝑣𝑥3(𝑡 − 𝑡2), 𝑡2 ≤ 𝑡 < 𝑡3

𝑥3, 𝑡 ≥ 𝑡3

(𝟕. 𝟏)

Com:

{

𝑑1 = |(𝑥1, 𝑦1, 𝑧1)|𝑑2 = |(𝑥1, 𝑦1, 𝑧1) − (𝑥2, 𝑦2, 𝑧2)|𝑑3 = |(𝑥2, 𝑦2, 𝑧2) − (𝑥3, 𝑦3, 𝑧3)|

{

𝑡1 =

𝑑1𝑣

𝑡2 = 𝑡1 +𝑑2𝑣

𝑡3 = 𝑡2 +𝑑3𝑣

{

𝑣𝑥1 = 𝑣

𝑥1𝑑1

𝑣𝑥2 = 𝑣𝑥2 − 𝑥1𝑑2

𝑣𝑥3 = 𝑣𝑥3 − 𝑥2𝑑3

(𝟕. 𝟐)

0

2

4

6

8

10

12

0 5 10 15 20 25 30 35

x (

m)

t (s)

xcomando

0

2

4

6

8

10

0 5 10 15 20 25 30 35

y (

m)

t (s)

ycomando

0

2

4

6

8

10

0 5 10 15 20 25 30 35

z (m

)

t (s)

zcomando

Em seguida, há o Controlador de Posição. Suas equações já foram mostradas no capítulo 6.1, só

vale ressaltar aqui alguns parâmetros adotados: 𝑓𝑚í𝑛 = 0,15𝑃, 𝑓𝑚á𝑥 = 2𝑃 e 𝛾𝑚á𝑥 = 30°. Na

verdade, 𝑓𝑚á𝑥 e 𝛾𝑚á𝑥 são relacionados:

cos(𝛾𝑚á𝑥) =𝑃

𝑓𝑚á𝑥 (𝟕. 𝟑)

Assim, dado 𝑓𝑚á𝑥 = 2𝑃, obtém-se 𝛾𝑚á𝑥 = 60°. Não há prejuízo em se fixar 𝛾𝑚á𝑥 abaixo

de 60°, somente acima.

As equações para o Controlador de Atitude também já foram exibidas, conforme o capítulo

6.2.

Em um modelo puramente teórico, as saídas dos controladores (𝑓, 𝜏𝑥, 𝜏𝑦 e 𝜏𝑧) poderiam

ser diretamente ligadas ao bloco dinamica. No entanto, na prática, o comando que aciona cada

motor é um throttle. Assim, são realizadas algumas conversões, conforme a Fig. 7.2.

A função torques2forca inverte o sistema definido pelas 𝐸𝑞. (3.5) e 𝐸𝑞. (6.2) para se

obter:

{

𝑓1 =

1

4(𝑓 + 𝜏𝑥

√2

𝑙− 𝜏𝑦

√2

𝑙+ 𝜏𝑥

1

𝑘)

𝑓2 =1

4(𝑓 − 𝜏𝑥

√2

𝑙− 𝜏𝑦

√2

𝑙− 𝜏𝑥

1

𝑘)

𝑓3 =1

4(𝑓 − 𝜏𝑥

√2

𝑙+ 𝜏𝑦

√2

𝑙+ 𝜏𝑥

1

𝑘)

𝑓4 =1

4(𝑓 + 𝜏𝑥

√2

𝑙+ 𝜏𝑦

√2

𝑙− 𝜏𝑥

1

𝑘)

(𝟕. 𝟒)

Essas forças de comando, por sua vez, são convertidas para as correspondentes velocidades

angulares do rotores, que por sua vez são convertidas para um throttle, que por fim aciona o motor

e gera as forças de comando:

𝑓 = 𝑘𝑎𝜔𝑚2 (𝟕. 𝟓)

𝑡ℎ𝑟𝑜𝑡𝑡𝑙𝑒 = 𝜔𝑚 (𝟕. 𝟔)

A 𝐸𝑞. (7.5) vem das leis aerodinâmicas e a 𝐸𝑞. (7.6) é uma mera simplificações, com o

objetivo somente de ilustrar a conversão para throttle.

Mais um motivo que torna importante explicitar as etapas intermediárias da Fig. 7.2 é o

fato de haver controladores entre os rotores e o computador de bordo. Basicamente, quando é

enviado um comando de voltagem para um motor, é um sistema de controle que faz a voltagem do

motor se aproximar do comando de voltagem proveniente do computador de bordo. Apesar deste

estudo não ter sido abordado neste trabalho, é importante mencionar que ele existe.

Por fim, para que o movimento obtido não esteja somente na forma de gráficos, adicionou-

se um bloco Visualização Simulink 3D Animation que anima em 3D o quadricóptero e sua

trajetória. Antes dele há o bloco m2D que converte os parâmetro MRP para a matriz de atitude

utilizando a 𝐸𝑞. (3.3).

8 Resultados e Discussão

A malha fechada foi simulada, sempre com condições iniciais nulas, como se o

quadricóptero estivesse inicialmente em repouso no chão. Os ganhos iniciais dos controladores

(𝑘1, 𝑘2, 𝑘3, 𝑘4, 𝑘𝑝 𝑒 𝑘𝑑) escolhidos foram arbitrários. Para o Controlador de Atitude, é prática

comum utilizar-se 𝑘𝑑 = 10%𝑘𝑝 para um primeiro teste, já que o ganho derivativo é mais

“agressivo”.

O procedimento realizado foi ajustar os ganhos buscando-se atingir os seguintes objetivos:

tempo de subida não muito longo, trajetória suave e overshoot pequeno.

Como, após o primeiro teste, os resultados não satisfazeram os objetivos, como era de se

esperar, os ganhos foram alterados. Um ganho era alterado por simulação para ser possível analisar

seu efeito isoladamente.

Sabe-se da teoria que quanto maior o ganho proporcional, mais rápida é a resposta, mas

também maior o overshoot. O ganho proporcional não pode ser muito pequeno porque isso tornaria

o sistema lento e pouco responsivo a distúrbios. Quanto maior o ganho derivativo, menor o

overshoot, mas ele pode aumentar as oscilações caso seja grande. Ele também não pode ser muito

pequeno, se não o overshoot causado pelo ganho proporcional não seria aliviado. Todas essas

relações foram de fato observadas durante as simulações, tendo sido utilizadas para ajustar os

ganhos de forma a atingir os objetivos de desempenho descritos.

Seguem os ganhos que resultaram em boas respostas do sistema:

𝑘1 = 1,5

𝑘2 = 2

𝑘3 = 1,05

𝑘4 = 1,5

𝑘𝑑 = 0,1

𝑘𝑝 = 1 (𝟖. 𝟏)

Para a simulação final, o quadricóptero deveria sair do ponto (0,0,0), ir para (7,8,9),

depois para (10,4,5) e retornar para o ponto de início (0,0,0). Preferiu-se realizar três trajetos logo

em um teste para fazer somente um teste, ao invés de realizar três testes diferentes.

Na Fig. 8.1, segue o rendimento em 3D da simulação, e, nas Fig. 8.2 a 8.13, os gráficos

referentes a esse trajeto.

Fig. 8.1 – Imagem em perspectiva da trajetória final

Fig. 8.2 – Final | 𝒙𝑪𝑴 Fig. 8.3 – Final | 𝒚𝑪𝑴 Fig. 8.4 – Final | 𝒛𝑪𝑴

Fig. 8.5 – Final | �̇�𝑪𝑴 Fig. 8.6 – Final | �̇�𝑪𝑴 Fig. 8.7 – Final | �̇�𝑪𝑴

Fig. 8.8 – Final | 𝜙 Fig. 8.9 – Final | 𝜽 Fig. 8.10 – Final | 𝝍

0

2

4

6

8

10

12

0 6 14 21 28 36

xC

M (

m)

t (s)

xCM

0

2

4

6

8

10

0 6 14 21 28 36

yC

M (

m)

t (s)

yCM

0

2

4

6

8

10

0 6 14 21 28 36

z CM

(m

)

t (s)

zCM

-1,5

-1,0

-0,5

0,0

0,5

1,0

0 6 14 21 28 36

vxC

M (

m/s

)

t (s)

vxCM

-1,0

-0,5

0,0

0,5

1,0

0 6 14 21 28 36

vyC

M (

m/s

)

t (s)

vyCM

-1,0

-0,5

0,0

0,5

1,0

0 6 14 21 28 36v

zCM

(m

/s)

t (s)

vzCM

-4

-2

0

2

4

6

0 6 14 21 28 36

Ro

ll (ᵒ)

t (s)

Roll

-6

-4

-2

0

2

4

0 6 14 21 28 36

Pit

ch (ᵒ)

t (s)

Pitch

-1,0

-0,5

0,0

0,5

1,0

0 6 14 21 28 36Ya

w (ᵒ)

t (s)

Yaw

Fig. 8.11 – Final | 𝝎𝒙 Fig. 8.12 – Final | 𝝎𝒚 Fig. 8.13 – Final | 𝝎𝒛

Nota-se, nas Fig. 8.5 a 8.7, que as componentes da velocidade são menores que 𝑣 = 1 𝑚/𝑠

e são praticamente constantes entre um ponto e outro a ser visitado, como foi fixado pela entrada

em rampa. Ainda, os ângulos de roll e pitch não ultrapassaram o limite 𝛾𝑚á𝑥 = 30°.

Comparando-se as Fig. 8.2, 8.3 e 8.4 às Fig.7.3, 7.4 e 7.5 vê-se que a trajetória obtida

condiz com a desejada. Ampliando-se os últimos instantes das Fig. 8.2 à 8.4, vê-se que o overshoot

é pequeno e a oscilação é mínima:

Fig. 8.14 – Ampliação da Fig. 8.2 Fig. 8.15 – Ampliação da Fig. 8.3 Fig. 8.16 – Ampliação da Fig. 8.4

Portanto, os resultados satisfizeram todas as imposições e objetivos feitos à simulação.

-0,1

-0,1

0,0

0,1

0,1

0,2

0 6 14 21 28 36ωx (

rad

/s)

t (s)

ωx

-0,2

-0,1

-0,1

0,0

0,1

0,1

0,2

0 6 14 21 28 36

ωy (

rad

/s)

t (s)

ωy

-1,0

-0,5

0,0

0,5

1,0

0 6 14 21 28 36

ωz

(ra

d/s

)

t (s)

ωz

-0,10

-0,05

0,00

0,05

0,10

0,15

0,20

0,25

34 36 38

xC

M (

m)

t (s)

xCM

-0,04

-0,02

0,00

0,02

0,04

0,06

0,08

0,10

34 36 38

yC

M (

m)

t (s)

yCM

-0,03

0,00

0,03

0,06

0,09

0,12

0,15

34 36 38

z CM

(m

)

t (s)

zCM

9 Conclusões e Sugestões para Trabalhos Posteriores

Quanto à parte das equações dinâmicas, houve dificuldades para determinar a melhor

parametrização para atitude e a equação cinemática para esta parametrização. Essas equações não

são triviais, tornando importante a escolha apropriada de sistemas coordenados e de parametrização

de atitude.

Quanto à parte do controle, foi possível perceber que ele possui particularidades

dependendo do sistema físico ao qual se aplica – neste estudo, por exemplo, há a questão da

saturação. O ajuste dos ganhos depende do desempenho que se deseja obter e da matriz de inércia

do corpo. No geral, o controle é fundamental para manter o sistema estável no caso da ação de

forças externas imprevisíveis e até mesmo no caso de simplificação das equações dinâmicas.

O mais importante foi que se cumpriu o objetivo de aprender sobre a dinâmica e o controle

de quadricóptero e de realizar uma simulação deste modelo com sucesso.

Para próximos trabalhos, seria interessante atribuir ao quadricóptero a missão de

transportar pequenos objetos. Seria necessário, neste caso, modificar o modelo físico, adicionando-

se uma garra robótica, por exemplo, e o objeto à matriz de inércia. Com a alteração do modelo

físico, os ganhos dos controladores deveriam ser reajustados.

Caso se desejasse construir o quadricóptero, os momentos de inércia seriam calculados em

laboratório por equipamentos próprios para identificar os eixos principais de inércia do corpo. Toda

a eletrônica do computador de bordo deveria ser implementada, especificando-se o funcionamento

dos blocos intermediários da Fig. 7.2.

Para medir o estado atual do corpo, seriam necessários odômetros e sensores. Estes podem

ser internos, tais como giroscópio, magnetômetros, sensores de pressão para altitude, ou externos,

tal como um conjunto de câmeras em uma sala. O primeiro teste prático seria em um globo que fixa

o centro de massa do corpo para verificar se a atitude comporta-se como esperado. Em um teste

final, o corpo teria seus seis graus de liberdade.

Referências

[1] SANTOS, Davi. Modelagem Dinâmica de um Quadricóptero.

[2] GYROFLY. Gyro200ED X4 - Manual do Usuário.

[3] SANTOS, Davi; SAOTOME, Osamu; CELA, Arben. Trajectory Control of Multirotor

Helicopters with Thrust Vector Constrains. June 2013.

[4] COSTA, Joana D. D. Projeto de um Controlador de Atitude de Baixo Custo para o Satélite

ITASAT-1. 2012.

Apêndice - Códigos do MATLAB

forces

function f = forces(t)

m = 1.2;

g = 9.81;

P = m*g;

% hovering

% f1 = P/4;

% f2 = P/4;

% f3 = P/4;

% f4 = P/4;

% movement along z axis (upwards)

% f1 = P/2;

% f2 = P/2;

% f3 = P/2;

% f4 = P/2;

% movement along z axis (downwards)

% f1 = P/8;

% f2 = P/8;

% f3 = P/8;

% f4 = P/8;

% movement along x axis (pitch=30 deg) m0=[0 1 0]tan(30/4)

% movement along y axis (roll=-30 deg) m0=[1 0 0]tan(-30/4)

% movement along rotor 1 (-30<pitch<0, 0<roll<30) m0=[-V2/2 V2/2 0]tan(30/4)

% f1 = P*sqrt(3)/6;

% f2 = P*sqrt(3)/6;

% f3 = P*sqrt(3)/6;

% f4 = P*sqrt(3)/6;

% yaw test (w_z > 0)

f1 = 3*P/8;

f2 = P/8;

f3 = 3*P/8;

f4 = P/8;

f = [f1 f2 f3 f4];

dinamica

function [sys,x0,str,ts,simStateCompliance]=dinamica(t,x,u,flag,xi,J,massa,g,l,k)

%tic

switch flag

%%%%%%%%%%%%%%%%%%

% Initialization %

%%%%%%%%%%%%%%%%%%

case 0

[sys,x0,str,ts,simStateCompliance] = mdlInitializeSizes(xi);

%%%%%%%%%%%%%%%

% Derivatives %

%%%%%%%%%%%%%%%

case 1

sys = mdlDerivatives(t,x,u,J,massa,g,l,k);

%%%%%%%%%%%%%%%%%%%%%%%%

% Update and Terminate %

%%%%%%%%%%%%%%%%%%%%%%%%

case {2,9}

sys = []; % do nothing

%%%%%%%%%%

% Output %

%%%%%%%%%%

case 3

sys = mdlOutputs(t,x,u);

otherwise

DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end

% end limintm

%

%=============================================================================

% mdlInitializeSizes

% Return the sizes, initial conditions, and sample times for the S-function.

%=============================================================================

%

function [sys,x0,str,ts,simStateCompliance] = mdlInitializeSizes(xi)

sizes = simsizes;

sizes.NumContStates = 12;

sizes.NumDiscStates = 0;

sizes.NumOutputs = 12;

sizes.NumInputs = 4;

sizes.DirFeedthrough = 0;

sizes.NumSampleTimes = 1;

sys = simsizes(sizes);

str = [];

x0 = xi;

ts = [0 0]; % sample time: [period, offset]

% specify that the simState for this s-function is same as the default

simStateCompliance = 'DefaultSimState';

% end mdlInitializeSizes

%

%=============================================================================

% mdlDerivatives

% Compute derivatives for continuous states.

%=============================================================================

%

function sys = mdlDerivatives(t,x,u,J,massa,g,l,k)

%% Declaração de Variáveis

% Vetor de estados [x y z vx vy vz m1 m2 m3 wx wy wz]

v = x(4:6);% velocidades lineares

m = x(7:9);% MRP

w = x(10:12);% velocidades angulares

f = u(1:4);

%% rotation_dynamics

% cálculo dos torques

Tx_c = l*sqrt(2)/2*(+f(1)-f(2)-f(3)+f(4));

Ty_c = l*sqrt(2)/2*(-f(1)-f(2)+f(3)+f(4));

Tz_c = k*(+f(1)-f(2)+f(3)-f(4));

% derivadas

wx_dot = ((J(2)-J(3))/J(1))*w(2)*w(3)+Tx_c/J(1);

wy_dot = ((J(3)-J(1))/J(2))*w(1)*w(3)+Ty_c/J(2);

wz_dot = ((J(1)-J(2))/J(3))*w(1)*w(2)+Tz_c/J(3);

%% attitude_kinematics

mRot = [0 -m(3) m(2); m(3) 0 -m(1); -m(2) m(1) 0];

m_dot = .25*((1-m'*m)*eye(3)+2*mRot+2*m*(m'))*w;

%% translation_dynamics

% rotation matrix (normalized)

D = eye(3) + (1/(1+m'*m)^2)*(8*mRot*mRot-4*(1-m'*m)*mRot);

f_tot = f(1)+f(2)+f(3)+f(4);

v_dot = (1/massa)*([0;0;-massa*g]+D'*[0;0;f_tot]);

sys = [v; v_dot; m_dot; wx_dot; wy_dot; wz_dot];

%end

% end mdlDerivatives

%

%=============================================================================

% mdlOutputs

% Return the output vector for the S-function

%=============================================================================

%

function sys = mdlOutputs(t,x,u)

sys = x(1:12);

%toc

% end mdlOutputs

m2angulosEuler

function angulos = m2angulosEuler(m)

% matriz de rotação mRot = [0 -m(3) m(2); m(3) 0 -m(1); -m(2) m(1) 0]; D = eye(3) + (1/(1+m'*m)^2)*(8*mRot*mRot-4*(1-m'*m)*mRot);

th = -asind(-D(3,1)); % pitch phi = -atan2d(D(3,2),D(3,3)); % roll psi = -atan2d(D(2,1),D(1,1)); % yaw %end

angulos = [phi; th; psi];

rampas

function r_d = rampas(t)

% waypoints

x1 = 7; x2 = 10; x3 = 0;

y1 = 8; y2 = 4; y3 = 0;

z1 = 9; z2 = 5; z3 = 0;

% distâncias

d1 = norm([x1-x3 y1-y3 z1-z3]);

d2 = norm([x1-x2 y1-y2 z1-z2]);

d3 = norm([x2-x3 y2-y3 z2-z3]);

% velocidades

v = 1; % m/s

vx1 = v*x1/d1;

vy1 = v*y1/d1;

vz1 = v*z1/d1;

vx2 = v*(x2-x1)/d2;

vy2 = v*(y2-y1)/d2;

vz2 = v*(z2-z1)/d2;

vx3 = v*(x3-x2)/d3;

vy3 = v*(y3-y2)/d3;

vz3 = v*(z3-z2)/d3;

xVec = zeros(381,1);

yVec = zeros(381,1);

zVec = zeros(381,1);

% tempos

t1 = d1/v;

t2 = t1 + d2/v;

t3 = t2 + d3/v;

if t < 0.2

x = 0.1;

y = 0.1;

z = 0.1;

else

if (0.2 <= t)&& (t < t1)

x = vx1*t;

y = vy1*t;

z = vz1*t;

else

if (t1 <= t) && (t < t2)

x = x1 + vx2*(t-t1);

y = y1 + vy2*(t-t1);

z = z1 + vz2*(t-t1);

else

if t <= t3

x = x2 + vx3*(t-t2);

y = y2 + vy3*(t-t2);

z = z2 + vz3*(t-t2);

else

x = x3;

y = y3;

z = z3;

end

end

end

end

r_d = [x y z];

rampasTeste

function rampasTeste

% waypoints

x1 = 7; x2 = 10; x3 = 0;

y1 = 8; y2 = 4; y3 = 0;

z1 = 9; z2 = 5; z3 = 0;

% distâncias

d1 = norm([x1-x3 y1-y3 z1-z3]);

d2 = norm([x1-x2 y1-y2 z1-z2]);

d3 = norm([x2-x3 y2-y3 z2-z3]);

% velocidades

v = 1; % m/s

vx1 = v*x1/d1;

vy1 = v*y1/d1;

vz1 = v*z1/d1;

vx2 = v*(x2-x1)/d2;

vy2 = v*(y2-y1)/d2;

vz2 = v*(z2-z1)/d2;

vx3 = v*(x3-x2)/d3;

vy3 = v*(y3-y2)/d3;

vz3 = v*(z3-z2)/d3;

xVec = zeros(381,1);

yVec = zeros(381,1);

zVec = zeros(381,1);

% tempos

t1 = d1/v;

t2 = t1 + d2/v;

t3 = t2 + d3/v;

i = 1;

for t = 0:0.1:38

if t < 0.2

x(i) = 0.1;

y(i) = 0.1;

z(i) = 0.1;

else

if (0.2 <= t)&& (t < t1)

xVec(i) = vx1*t;

yVec(i) = vy1*t;

zVec(i) = vz1*t;

else

if (t1 <= t) && (t < t2)

xVec(i) = x1 + vx2*(t-t1);

yVec(i) = y1 + vy2*(t-t1);

zVec(i) = z1 + vz2*(t-t1);

else

if t <= t3

xVec(i) = x2 + vx3*(t-t2);

yVec(i) = y2 + vy3*(t-t2);

zVec(i) = z2 + vz3*(t-t2);

else

xVec(i) = x3;

yVec(i) = y3;

zVec(i) = z3;

end

end

end

end

i = i+1;

end

tempo = 0:0.1:38;

plot(tempo,xVec);

% plot(tempo,yVec);

% plot(tempo,zVec);

controladorPosicao

function saidas = controladorPosicao(entrada)

% entradas

r_d = entrada(1:3); % posicao desejada

r = entrada(4:6);

v = entrada(7:9);

m = entrada(10:12);

% constantes

k1 = 1.5;

k2 = 2;

k3 = 1.05;

k4 = 1.5;

gama_max = 30*pi/180;

massa = 1.2;

g = 9.81;

P = massa*g;

f_min = 0.15*P;

f_max = 2*P;

% calculo de n3 e gama

mRot = [0 -m(3) m(2); m(3) 0 -m(1); -m(2) m(1) 0];

D = eye(3) + (1/(1+m'*m)^2)*(8*mRot*mRot-4*(1-m'*m)*mRot);

n3 = D(3,3);

% calculo de f

f_esperado = (massa/n3)*(g-k1*(r(3)-r_d(3))-k2*v(3));

if (f_esperado < f_min)

f = f_min;

else

if (f_esperado > f_max)

f = f_max;

else

f = f_esperado;

end

end

% calculo de n12

n12_esperado = (-massa/f)*(k3*(r(1:2)-r_d(1:2)) + k4*v(1:2));

if norm(n12_esperado) <= sin(gama_max)

n12 = n12_esperado;

else

n12 = sin(gama_max)*(n12_esperado/norm(n12_esperado));

end

% n12 para angulo e eixo principal de Euler

dir = cross([0; 0; 1],[n12; 0]);

e = dir/norm(dir);

gama_d = asin(norm(n12));

% angulo e eixo principal de Euler para MRP

m_d = tan(gama_d/4)*e;

saidas = [f; m_d];

controladorAtitude

function T_c = controladorAtitude(entradas)

% entradas

m_d = entradas(1:3);

m = entradas(4:6);

w = entradas(7:9);

% ganhos

Kp = 1;

Kd = 0.10*Kp;

% lei de controle

T_c = Kp*(m_d-m) - Kd*w;

torques2forca

function F = torques2forca(T)

l = 0.277; % braco

k = 0.01;

f = T(1);

Tx_c = T(2);

Ty_c = T(3);

Tz_c = T(4);

a = sqrt(2)/l;

f1 = .25*(f + Tx_c*a - Ty_c*a + Tz_c/k);

f2 = .25*(f - Tx_c*a - Ty_c*a - Tz_c/k);

f3 = .25*(f - Tx_c*a + Ty_c*a + Tz_c/k);

f4 = .25*(f + Tx_c*a + Ty_c*a - Tz_c/k);

F = [f1 f2 f3 f4];

f2w

function W = f2w(F)

k = 1;

W1 = sqrt(F(1)/k);

W2 = sqrt(F(2)/k);

W3 = sqrt(F(3)/k);

W4 = sqrt(F(4)/k);

W = [W1;W2;W3;W4];

w2th

function th = w2th(W)

th = W;

th2f

function f = th2f(th)

k = 1;

f1 = k*th(1)^2;

f2 = k*th(2)^2;

f3 = k*th(3)^2;

f4 = k*th(4)^2;

f = [f1;f2;f3;f4];

m2D

function D = m2D(m)

mRot = [0 -m(3) m(2); m(3) 0 -m(1); -m(2) m(1) 0];

% rotation matrix (normalized)

D = eye(3) + (1/(1+m'*m)^2)*(8*mRot*mRot-4*(1-m'*m)*mRot);

FOLHA DE REGISTRO DO DOCUMENTO 1. CLASSIFICAÇÃO/TIPO

TC

2. DATA

20 de novembro de 2014

3. REGISTRO N°

CTA/ITA/TC-077/2014

4. N° DE PÁGINAS

57 5. TÍTULO E SUBTÍTULO:

Simulação de um sistema de controle de um quadricóptero. 6. AUTOR(ES):

Tarik Carvalho Bussab 7. INSTITUIÇÃO(ÕES)/ÓRGÃO(S) INTERNO(S)/DIVISÃO(ÕES):

Instituto Tecnológico de Aeronáutica – ITA 8. PALAVRAS-CHAVE SUGERIDAS PELO AUTOR:

Quadricóptero, Simulação, Controlador de Atitude, Controlador de Posição, Dinâmica, Modelagem. 9.PALAVRAS-CHAVE RESULTANTES DE INDEXAÇÃO:

Controle de atitude; Simulação computadorizada; Modelagem; Dinâmica; Aeronaves; Engenharia

aeronáutica. 10.

APRESENTAÇÃO: X Nacional Internacional

ITA, São José dos Campos. Curso de Graduação em Engenharia Mecânica-Aeronáutica. Orientador: Prof.

Davi Antônio dos Santos. Publicado em 2014.

11. RESUMO:

Este trabalho se dedicará à modelagem, simulação e controle de atitude e posição de um robô aéreo do

tipo quadricóptero. É bem conhecido que o veículo em questão apresenta uma dinâmica subatuada, pois

contém seis graus de libertade (três de rotação e três de translação), enquanto está sujeito a apenas quatro

variáveis de atuação independentes (três componentes de torque e um componente de força perpendicular

ao plano dos rotores). Devido a essa característica, o quadricóptero não poderia ser controlado

independentemente em seus seis graus de liberdade. Sem embargo, isso não consiste numa dificuldade,

desde que o interesse central é normalmente controlar independentemente apenas quatro graus de

liberdade: os três componentes de posição e o ângulo de heading. O sistema que realiza tal controle é

composto de duas malhas aninhadas, em que a interna se ocupa do controle de atitude e a externa é

responsável por controlar a posição. A modelagem dinâmica será realizada mediante a formulação

Newton-Euler. Na malha interna (de atitude) far-se-á uso de três controladores proporcionais-derivativos

(PD) desacoplados. O controle da malha externa (de posição) será realizado pela lei de controle proposta

em [3], que consiste num controlador PD saturado em combinação com uma linearização por

realimentação de estados. O sistema descrito será simulado no ambiente MATLAB/Simulink. A

simulação consiste em o robô ser comandado a passar por dois pontos no espaço e depois retornar à sua

posição inicial.

12. GRAU DE SIGILO:

(X ) OSTENSIVO ( ) RESERVADO ( ) CONFIDENCIAL ( ) SECRETO